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ABSTRACT 

The local theory of the critical lines of 2D and 3D Gaussian fields that underline 
the cosmic structures is presented. In the context of cosmological matter distribution 
the subset of critical lines of the 3D density field serves to delineate the skeleton 
of the observed filamentary structure at large scales. A stiff approximation used to 
quantitatively describe the filamentary skeleton shows that the flux of the skeleton 
lines is related to the average Gaussian curvature of the one D sections of the field, 
much in the same way as the density of the peaks. The distribution of the length of the 
critical lines with threshold is analyzed in detail, while the extended descriptors of the 
skeleton - its curvature and its singular points, are introduced and briefiy described. 
Theoretical predictions are compared to measurements of the skeleton in realizations of 
Gaussian random fields in 2D and 3D. It is found that the stiff approximation predicts 
accurately the shape of the differential length, allows for analytical insight, and explicit 
closed form solutions. Finally, it provides a simple classification of the singular points 
of the critical lines: i) critical points; ii) bifurcation points; iii) slopping plateaux. 



1 INTRODUCTION 



The concept of random fields is central to cosmology. Random fields both provide initial conditions for the evolution of the 
matter distribution in the Universe, and represent how the observed signals manifest themselves in 3D, (e.g., in the galaxy or 
matter density inhomogeneities that form the Large Scale Structure (LSS)), or on the 2D sky (e.g. for the Cosmic Microwave 
Background (CMB) temperature and polarization, the convergence or shear in weak lensing maps). In the modern cosmological 
theories where initial seeds for inhomogeneities observed as cosmic structures have quantum origin, the fields of initial density 
fluctuations (and velocities) are Gaussian. Subsequent evolution retains Gaussianity for the observables that evolve linearly 
(CMB, very Large Scale Structure) while developing non-Gaussian signature if non-linear effects are involved (e.g. lensing 
and LSS at smaller scales). 

While comparing the observational data to cosmological theory, in particular in order to estimate parameters of cosmo- 
logical models, the emphasis is traditionally placed on the statistical descriptors of the random fields. For Gaussian fields 
the two-point correlation function or the power spectrum provide full statistical information, while non-Gaussian properties 
may be reflected in multi-point correlations. The understanding and the description of the morphology of structures in our 
Universe, on the other hand, calls for the studies of the geometry and topology of random fields. This subject has an extensive 
history from the early description of the one-dimensional radio signal time-streams in 1940's, to the study of the 2D ocean 
wave patterns in 1960's ( Longuet-Higgins|T957|) to 3D d imensional fields ( Adler|1981 1 that found the most fruitful application 
in cosmology ( Arnol'd et al.|1981 Bardeen et al.|1986 |. The most prominent geometrical objects in a typical realization of a 
random field are rear events - regions of unusually high or low values of the field. The rare events are usually related to the 
most spectacular observed objects - clusters of galaxies at low z, large protogalaxies at high-2; or extensive voids. They are 
associated with the neighbourhoods of extrema - maxima or minima - making studies of such critical points the first step in 
understanding typical geometry of a field ( |Kaiser|1984 , Bardeen et al.|[l986[ [Regos fc Szalay||1995[ [Scannapieco et al.|[2006| ). 
The behaviour of the field in the neighbourhood of a rare peak is highly correlated with the peak properties, which allows to 
describe not only extrema but the extended peak-patch region ( Bond fc Myers||1996a I as a point process that involves the 
field and its successive derivatives. Including the shear flow into consideration gives a compelling application of the geometry 
of rare events to the description of cluster formation through the peak-patch collapse ( Bond fc Myers|[l996b I. 
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The rare events reflect the organization of the fleld around them and by and large determine the way the high (low) field 
regions are interconnected by the bridges of enhanced field values. In application to cosmology, the "Cosmic Web" picture 
emerges, which relates the observed clusters of galaxies, and filaments that link them, to the geometrical properties of the 
initial density field that are enhanced but not yet destroyed by the still mildly non-linear evolution on supercluster scales 
( Bond et al.||1996 1 . The study of the connectivity of filamentary structures reveal the role of the remaining type of critical 
points, the saddle extrema, in establishing, in particular, the percolation properties of the Web ( Colombi et al.|2000 |. The next 
step naturally involves describing the statistical properties of these filamentary structures ( Pogosyan et al.|1998 Schmalzing 



et al.|1999 1 and developing techniques for mapping the filaments in the simulation and data. Novikov et al. (j2006) presented 



a 2D algorithm to trace the filaments of a density field while introducing the skeleton as the set of locally defined critical lines 
emanating from the critical points. Sousbie et al. (20081 (hereafter SPCNP) extended the local theory and algorithm to three 
dimensions and provided the foundation for this work while introducing the "stifi^" approximation. Recently, [Sousbie et al.| 
( 2008 1 presented an algorithm to map out a fully connected version of the skeleton that is defined according to the global 
properties as the lines of intersections of the patches (see also Aragon-Calvo et al. (20071, Platen et al. (20071 for alternative 
algorithms). This approach connects the study of the filamentary structure to the geometrical and topological aspects of the 
theory of gradient fiows ( Jost|2008 l and returns the focus to the notions of peak and void patches. 

This paper presents a consistent local theory for the cosmic skeleton, while focusing on the stiff approximation to compute 
the differential length of the skeleton as a function of the contrast and modulus of the gradient of the field. It allows us to 
define precisely how the properties of the skeleton depend analytically on the underlying spectral parameters, and understand 
what type of line prevails where. The crucial advantage of the local approach to the critical lines is that it allows to cast 
the statistical treatment of the linear objects as a point process that involves the field and its derivatives, which allows for 
analytical insight, and explicit closed form solutions. Our purpose is to construct the theory of critical lines of a given field 
corresponding to an intermediate representation of the field, which is more extended than the knowledge of the critical points. 

The organization of the paper is the following. Section [2] classifies the various critical lines in 2D and 3D, connects the 
average length in a unit volume to the flux of the skeleton lines and, within the stiff approximation, to the average Gaussian 
curvature of the field in transverse sections. It also discusses the meaning of this approximation. Section |3] calculates the 
differential length of all sets of critical lines in 2D, while Section |4] investigates the corresponding 3D set of critical lines. 
More generally, the expression for the differential length of the N dimensional skeleton is sketched in Appendix]^ Section [5] 
introduces the extended descriptors of the skeleton, Section [5.3| describes their singular points, while Section [5] provides the 
discussion and the summary. Appendix|D]gives the general method for obtaining in close form the joint distribution of the field 
and any combination of its derivative tensors in arbitrary dimensions. In particular, it exhibits all the statistical invariants 
and their dependence on the spectral parameters. 



2 THE CRITICAL LINES AND THE SKELETON OF A 3D RANDOM FIELD 
2.1 Local definition and classification 

The subject of our investigation is a random field, p(r), that in a cosmological setting describes, for example, the density 
of the matter in the Universe, or the projected distribution of Cosmic Microwave light on the celestial sphere. Our focus is 
on the geometrical properties of the critical lines, that connect extrema of the fleld mapping out the fllamentary ridges and 
valleys of the field. SPCNP have introduced the definition of the local eritieal lines as the set of points where the gradient of 
the density, Vp, is an eigenvector of its Hessian matrix, Ti, Ti. ■ Vp = AVp i.e, the gradient and one of the principal curvature 
axes are coUinear. Formally, this can be specified by a set of equations 

S= (Vp-H)-e-(Vp) = 0, (1) 

where e is the fully antisymmetric (Levi-Civita) tensor of rank TV. In general S is an antisymmetric N — 2 tensor. 

In 3D, the function S is vector-valued, S" — e^'^\Vmp)'H"^k('^ip)- However, zeroes of S determine a set of lines 

rather than isolated points. Let us consider the behaviour of S function in the neighbourhood of a point r = that satisfies 
criticality condition S(0) = 0: 

S((5r) RiO-f^(VfcS)5r'= . (2) 
fe 

In our case, under the condition S" = 0, the matrix VfeS" by definition possesses the left null- vector, furnished by the 
density gradient, Vip (Vfc5") — 0; hence the gradients VS" are not linearly independent. Consequently, there is a non- 
trivial solution for the nght null-vector (VfcS") Sr'° = 0, which determines the local direction of the line along which the 
criticality condition is maintained, S(5r) — 0. The critical lines intersect where Vfe5" admits more than one independent right 
null- vector. 

When we take the eigenvalues of the Hessian to be sorted, Ai ^ A2 ^ A3, the gradient of the field at the critical line may 
be found aligned with the first, second, or third eigenvector. This gives rise to the classification of the critical lines based on the 
choice of the eigenvector aligned with the gradient, that becomes more fine grained when the magnitudes of the eigenvalues 
are taken into account. Namely, we distinguish primary critical lines, which correspond to Vp being aligned with the direction 
in which the field is the least curved, i.e where the eigenvalue is the smallest in magnitude, and secondary critical lines at 
which Vp is aligned with the eigenvalues of larger magnitude. The primary type consists of 
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Type 




Alignment 




Condition 




Primary Skeleton; 


Ti. 


• Vp = AiVp 


Ai - 


- A2 ^ 




Inter-skeleton: 


n 


. Vp = A2VP 


Ai - 


1- A2 > and A3 - 


h A2 < 


Anti-skeleton: 


n 


■ Vp = AgVp 


A3- 


1- A2 ^ 




Secondary 


n 


• Vp = A2VP 


Ai- 


1- A2 ^ 






n 


■ Vp = AgVp 


Ai- 


1- A2 ^ 






n 


. Vp = AiVp 


Ai- 


1- A2 > and A3 - 


h A2 < 




H 


• Vp = AgVp 


Ai- 


1- A2 > and A3 - 


h A2 < 




n 


. Vp = AiVp 


A3 - 


1- A2 ^ 






n 


• Vp = A2VP 


A3- 


1- A2 ^ 





Table 1. the classification of the critical lines in 3D. 



(i) The skeleton, that has the gradient in the Ai direction and is hmited to the region |Ai| ^ IA2I, which translates to the 
condition Ai -|- A2 ^ 0. The skeleton has always eigenvalues in the directions transverse to Vp negative, X3 ^ X2 ^ and 
corresponds to the filamentary ridges spreading from the maxima in the direction of the slowest descent. 

(ii) The anti-skeleton, that has the gradient in the A3 direction and is restricted to the region IA3I ^ IA2I, i.e A3 + A2 5S 0. In 
the directions transverse to Vp the anti-skeleton has always positive curvature Ai A2 0. It corresponds to the filamentary 
valleys spreading from the minima in the direction of the slowest ascent. Anti-skeleton can be viewed as a skeleton of the — p 
field. 

(iii) The intermediate skeleton along which the gradient is aligned with the middle eigen-direction of the Hessian where 
this direction is the shallowest IA2I < |Ai|, IA3I, i.c — Ai < A2 < —A3. This conditions is only possible in saddle- like regions 
where Ai > and A3 < 0. 

The formal classification of the critical lines is summarized in Table [T] 



2.2 The average flux (length per unit volume) of the critical lines 

As the average number density is the fundamental quantity that describes point events, e.g. extrema of a field, conversely 
the most important characterization of the critical lines or skeleton is their flux, i.e. the number of critical lines intersecting 
a given oriented surfac^ This flux is equivalent to the length of the lines per unit volume. Following NCD and SPCNP we 
shall preferentially use the latter terminology as it highlights that we deal with the first geometrical parameter, the length, 
of the lines. The subsequent parameters of these linear objects are the curvature, and, in 3D, the torsion. 

In this paper we consider p to be a homogeneous and isotropic Gaussian random field of zero mean, described by the 
power spectrum P{k). In the statistical description of the skeleton of the field p, several linear scales are involved 

-Ro = — , R, = — , it — — , R = — , (3) 

cri (T2 ""3 o"4 



where 



9 -D/2 roo 

'^o = (p'), a? = ((Vp)2), al = {{^pf), al = {(y^pf), and generally 4 = P{k)k''-^ dk . 

(4) 

These scales are ordered Ro ^ R, ^ R ^ . . .. The first two have well-known meanings of typical separation between zero- 
crossing of the field -Ro and mean distance between extrema, ii, (Bardeen & al. 1986), and the third one, R is, by analogy, the 
typical distance between the infiection points. These three are the only ones that are involved in determination of the length 
of the critical lines. The higher order scale R appear only in computation of the curvature and the torsion (see Section [5|. 

Let us define a set of spectral parameters that depend on the shape of the underlying power spectrum. Out of these five 
scales four dimensionless ratios may be constructed that are intrinsic parameters of the theory 

7= — = , 7= — = , 7=^ = , and generally 7p,g = . (5) 

From the geometrical point of view 7 specifies how frequently one encounters a maximum between two zero crossings of the 
field, while 7 describes, on average, how many inflection points are between two extrema. From a statistical perspective, 7's 
are the cross-correlation coefficients between the field and its derivatives at the same point (see Appendix [D|. 

7=i^, 7^ (VP-VAP) ^ 



^ For M dimensional objects in N dimensional space, in general, one counts the average number of intersections between objects M and 
N-M dimensional surfaces, per unit N-M volume. From a statistical point of view this constitutes a point process that can be evaluated 
knowing the distribution of the field and some of its derivatives at one arbitrary point only. 
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For Gaussian fields, tliese parameters can be easily calculated from the power spectrum. All 7's range from zero to one. For ref- 
erence, for the power-law spectra with index n > —3, smoothed at small scales with a Gaussian window, 7 — ^(n -f 3)/(n + 5), 
7 — y''{n + 5)/{n + 7). Note that cosmologically relevant density power spectra have n > —3 and, thus, while 7 can attain 
low values, 7 are always close to unitj[^ . 

Let us introduce the dimensionless quantities for the field and its derivatives as well as for the functions 5™ and their 
gradients V5": 

aox = p, a-iXk = Vkp, a2Xki =\7kyip, a^Xkim = y N kp, (72(Tis' = 5', (T2(JiVs* = V5\ 0-40-1 VVs' = VV5', 
giving 

= ^ e'"''XmlXlXk , and VmS' = ^ e"'' {^~^ X^lmXlXk + [XnlXlmXk + X^lXkniXl]) . (8) 
klm kin 

Note the specific choice of scaling for VS which is convenient in view of the subsequent development of the so-called "stiff" 
approximation. SPCNP has shown that in terms of these dimensionless quantities, the cumulative length per unit volume of 
the total set of critical lines below the threshold i] is given by 

^-^^ J dx j d^Xkd^Xkld^°Xklni l^s' X \/s^\P{x,Xk,Xkl,Xklrn)SD (^s\Xk,Xkl,Xklm.)'^ So ^S^ {xk,Xkl,Xklni)'^ , (9) 

where a pair Vs\ Vs-* can be chosen arbitrarily as long as it is linearly independent. In this equation |Vs* x Vs-'l reflects the 
inverse characteristic area orthogonal to a critical line per one such line while the two (Jo-functions account for the critical line 
condition ([TJ. For the complete set of critical lines, there are no restriction to the region of integration. If one is interested in a 
particular type of the critical lines, the integration should be restricted to the regions consistent with Table[l] The differential 
length (per unit volume) is simply given by the derivative of equation ([9| with respect to rj: 

-r^ = ) J '^^^kd^^kld^^Xklm, |Vs' X '\/s^P{r!,Xk,Xkl,Xklm)SD [s\x k , X kl , X klm)^ 5d (^S^ {Xk , Xkl , Xklm)^ , (10) 

while the total length of critical lines is 

L = C{oo) = j dn^ . (11) 

Since for Gaussian field, the derivatives of even order are uncorrelated with the odd orders, the joint distribution function 
P{x, Xk, Xki, Xkim) entering equation Q is factorized as 

P{x,Xk,Xkl,Xklm) = Po{x,Xkl)Pl{Xk,Xklm)- (12) 

In Po, the only dependence on the power spectrum of the field is through the parameter 7 (c.f. equation ([5|) that describes the 
correlation between the field and its second derivatives. Similarly P\(xk,Xkim) only involves 7 which describes the correlation 
between the gradient of the field and its third derivatives. Therefore, dC/drj depends only on rj, R and 7. The integrated 
length, L may depend only on 7 and R since the marginalization of Po{ri,Xki) over rj eliminates the dependency over 7. 

2.3 The "stiff" filament approximation 

Let us look at the dependence of the length of the critical lines on characteristic scales of the field in more detail. The R~^ 



factor that appeared in equation ( 10 1 reflect our choice of dimensionless variables Q and is suggestive but not yet conclusive 



since jVs' x Vs-'| that includes third derivative terms, depends also on the other scale, R. Let us write formally 

Vs* X Vs^ = ^'^A{xk,Xki, Xkim) + 7"^B(a;fc, Xki,Xkim) + C{xkXki) ■ (13) 

If the third derivatives are important and the flrst term dominates, then the length scaling L oc ^~'^R~^ = R~^ would reflect 
the mean separation between the inflection points, R. Indeed, by deflnition the local skeleton is almost straight within a 
volume that has one inflection point ~ _R^. A straight segment through such volume has length ~ R, thus the expected length 
per unit volume is ~ 1/7?-^. But if the last term dominates statistically, the length per unit volume of the skeleton will scale as 
L oc RZ^ that can be interpreted that the critical lines are almost straight within a large volume volume ~ Ri containing one 
extremum. This is consistent with observation that the integral term does not depend on the third derivatives, thus inflection 
points play no role, and any dependence on 7 drops out. 

Which regime holds can be established by measuring the dependence of the critical lines length in the simulations as a 
function of smoothing length for different spectral indexes. For the power-law spectra with Gaussian smoothing at the radial 
scale a, in 3D, i?, — \f2a j \Jn 4- 5, while R = ^/2a/y/n + 7. The measurements in SPCNP found that L oc {n + 5.5)o~^ over 
the range of spectral indexes relevant to cosmology, which points at the subdominant nature played by the third derivatives. 
In the "stiff" approximation we omit the third derivative, effectively assuming that the Hessian can be treated as constant 

^ Cosmological density fields, therefore, have of order one inflection point per extremum, unUke, for, example, a mountain range, where 
one encounters many inflection points on a way from a mountain top to the bottom; see also Section 14.41 
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Figure 1. The neighbourhood of a local critical line (thick blue line). This is a zoomed section of the wide patch shown in Figure|2] Thin 
lines are isocontours of the field. Three sample points are investigated in detail. The signature, orientation and the magnitude of the local 
Hessian are represented by the golden shapes. Near the maximum on the right edge, the signature of the eigenvalues of the Hessian is 
(-1,-1), which is shown by ellipses oriented according to eigen-directions with longer semi-axis along the direction of the least curvature. 
At the leftmost point the eigenvalue signature is "saddle-like", (1,-1), which is represented by a pair of hyperbolae, also oriented with 
respect to eigen-directions. By definition, on the critical line the gradient of the field Vp, shown by red arrows, is aligned with one of the 
eigen-directions (i.e the axis of the ellipse or hyperbola in the graph). The light cyan arrows are the tangent vectors to the critical line 
oc e • VS, while stiff approximation to them would be parallel to the gradient. The direction of the critical line is close to the gradient 
when it follows the ridge near the maximum, but slides at an angle in the "saddle-like" region, before joining the saddle extremal point 
beyond the left edge of the plot (see Figure [2]l. Note that the gradient line that takes us to the same saddle as a segment of the global 
skeleton (dashed green) does not follow the ridge too closely in this instance. 



during the evaluation of Vs. This picture corresponds to a skeleton connecting extrema with relatively straight segments. In 



the "stiff" approximation, equation ( 10 \ becomes 



J d^Xk(f'Xkl \C{xk,Xlm)\ Po(ri,Xkl) Pl{Xk)SD (^s'{Xk,Xkl)'^ Sd (^S^Xk,Xkl)'^ . (14) 



The differential length is then only the function of 7 times L. 

The "stiff" approximation can be looked at from another perspective. By definition at a point on a local critical line, 
two of the characteristic directions defined for the field, namely, the direction of the gradient, Vp, and one chosen eigen 
direction of the Hessian, H, must coincide. But the direction of the critical line itself, given by V5" x VS-*, is not, in 
general, aligned with the gradient of the field. Local critical lines are not the gradient lines, and in this sense they differ from 



the skeleton lines defined globally as void-patch intersections (Sousbie et al. 20081. In the "stiff" approximation, however, 
(VmS') g ^ J^kin'^''"'' + XriiXkmXi] and [(Vs')_,y^ X (Vs^)^^.^J X Vp = 0, i.e. it is parallel to the gradient. 

Figure n] shows the details of the calculations for the high-resolution segment of the 2D field. Thus, the essence of the stiff 
approximation lies in the assumption that the mismatch between the critical lines and gradient directions is statistically small. 
As Figure |2] which contains an extended view of the same field, illustrates, this assumption holds particularly well for the 
primary critical lines which more closely correspond to the intuitive picture of sharp ridges and deep valleys. Indeed, at a 
primary line the gradient points to the least curved direction, i.e, in some sense in the direction in which the changes of the 
field properties are the slowest. Therefore one can expect that this is the direction in which the condition of criticality will 
be maintained, i.e which the critical line itself will follow. Figure [2] shows that the primary lines start to deviate from the 
gradient fiow mostly towards their end points when the curvature of the field along the line becomes comparable in magnitude 
to the transverse one. Secondary critical lines are much less certain to follow the gradient, sometimes exhibiting a "sliding" 
behaviour, on occasion almost orthogonal to the gradient, as a loop-like secondary line near the right saddle in Figure [2] 
exhibits. So the stiff approximation for the secondary lines should be taken with more caution, although we include them for 
completeness. 

The stiff approximation provides a framework to compute the total differential length of the critical lines and the local 
skeleton almost completely analytically. In the next two sections we will carry this calculation in two and three dimensions 
and argue that it can straightforwardly be extended in N dimensions (see Appendix [AJ . In what follows we shall omit in the 
derivation for brevity the 1/R* (in 2D) and (in 3D) factors, but keep in mind that all the length quantities below scale 

accordingly. In section [4.4] after the computational machinery is developed, we return to the role the third derivative may 
play in description of the critical lines. 
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Figure 2. An example of a generic patch of a 2D field. The underlying isocontours correspond to the density field. The thin gold lines 
arc the gradients lines of the field. The blue lines is the local set of critical lines, given by the solution of cS = {Ti. ■ Vp) X Vp = 0. 
Primary lines are shown in solid and the secondary lines are dashed. The green lines correspond to global critical lines: the skeleton and 
the anti skeleton, which delineate a special bundle of gradient lines | |Jost||2008j l at the intersection of a peak patch and a void patch. 
The primary local lines follow fairly well the gradient lines, noticeably near the extrema, where the "stiff" approximation holds best. 
In contrast, the approximation worsens for the secondary critical lines. The main distinction between the global and local skeletons is 
that the global one follows the everywhere smooth gradient line that uniquely connects a maximum to a saddle, at the cost of deviating 
from being exactly on the ridge (see how in the vicinity of the minimum at the bottom, the right green line does not follow the trough) 
The local skeleton tries to delineate the ridges as far from extrema as possible, but then the lines that follow this local procedure from 
different extrema do not meet and have to rather suddenly reconnect. A particularly striking example of this is the loop on the right 
hand side. A zoomed view of the area left to the top maximum is shown in Figure [l] 



Type Alignment Condition 



Primary 


Skeleton: 


n 


Vp = 


AiVp 


Ai - 


I-A2 ^0 




Anti-skeleton: 


H 


Vp = 


A2VP 


Ai- 


1- A2 > 


Secondary 




n 


Vp = 


A2VP 


Ai - 


I-A2 sJO 






H 


Vp = 


AiVp 


Ai- 


1- A2 > 



Table 2. the classification of the critical lines in 2D. 



3 CRITICAL LINES OF 2D FIELDS 

Even though the large scale structures of the universe are three dimensional, other important observed data sets involve 2D 
maps such as the cosmic microwave background or lensing convergence maps. Hence analyzing the local statistical properties 
of filaments in two dimensions is astrophysically well motivated. The 2D case is also a convenient starting point to introduce 
the details of the calculations that can be generalized to 3D and higher dimensions. 

The 2D case affords several simplifications over the 3D case. In 2D, 5 is a (pseudo) scalar function and its zero level, 
orthogonal to V S, determines the critical lines. The expression for the differential length simplifies to 

^ = y (fxk<fxkld'^Xklm \Vs\P{rj,Xk,Xkl,Xklm)SD {s{Xk,Xkl)) ■ (15) 

There are just four types of critical lines: two primary, the skeleton and the anti-skeleton, and two corresponding secondary 
ones. The classification of the 2D critical lines is summarized in Table [2] We shall focus on the most interesting primary lines 
in the main text, leaving the secondaries to the Appendix [B] In Figure|2]the critical lines of different types are shown for an 
example generic patch of a 2D field. 
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3.1 The differential length of the critical lines of 2D fields 

For 2D Gaussian fields, the calculation of the length of the critical lines can be carried almost completely analytically in the 
stiff approximation. 



3.1.1 Direct derivation in the field's frame 

Let us first proceed in the original coordinate frame. Defining 

S = XiX2 (Xli — X22) + X\2 {xl — xX) = IwXxX^ + 3^12 {x\ — X^) , 

the stiff approximation to Vs involves only up to second derivatives of the field 



I Vs|^ = [xl + xl) (w^ + XI2) ( + 4 (lu^ + Xr2 ) 



2a;iX2a;i2 -\- w [x-^ — X2) 
- 4 2 \ -u 

3j ~\~ 



and equation ( 15 1 becomes explicitly 
: exp 



dC 



16 



dv (27r)Vl - 7^ 



!L 
2 



JJJJJ dud'wdx\2dx\dx2 lVs| (5d(s) exp 



(u — 'yri)'^ ,2,2 2 2 

— ^-4™ - 4a;i2 - Xl - a;2 

2(1-72) 



(16) 



(17) 



(18) 



where the second derivatives are described using u — —{xn +2:22), w — {xn -~X22)/2 and X12. Let us integrate over X12 u sing 
the fc-function, which leads to a substitution X12 ^ {2xiX2w) / {x^ — x\) with the Jacobian \ l/{x\ ~x%)\. Then equation ( 18 1 
becomes 



dC 



16 



df) (27r)3^1 -7^ 



: exp 



' 2 



//// 



dudwdx\dx2 -, 



IVsl 



exp 



{'U.-l-nf . 2{xl+xlf 2 2 

- 4to ^TTT — X-i — X2 



2(1-7^) 



{xl-xlY 



where 



2 ^,2{xl + xlf ( , „ {xi + xi) 



IVsl = w 



Let us now substitut^E] 



(x\ X2 
(xl+xl) 



2x\ 2 



u + 2w 



{xi x\) 



{x\ X2 



to obtain 



dC 



16 



dv (27r)Vl 



exp 



noting that w' 

w{u + 2to)| 



(u)^ + XI2) , 



II II dud'wdxidx2 

JJJJ ^x\ + 



■ exp 



L '2(1-7^) 



. ~2 2 2 
iW — Xl — X2 



The integration over the first derivatives is now easily performed in the polar coordinates of the xi , X2 plane to give 



dC 



: exp 



2 



00 POO 



dudw \ w(u + 2w) \ Gxp 



oo-/ —00 



This is the final integral form which can be easily investigated in the u, w plane. 



4^^-4»^' 
2(1-72) 



(19) 



(20) 



(21) 



(22) 



(23) 



3.1.2 Derivation in the Hessian eigenframe 

To generalize the derivation to higher dimensions we note that one can just perform all the calculations in the Hessian 
eigenframe. We shall denote all quantities evaluatedin the eigenframe with tilde, e.g., a;ii(= Ai),3;22(= A2), ii, w, aii, i2. What 
must be taken into account is that, in general, these quantities are not Gaussian random variables (while the corresponding 
ones in the fixed frame are), since the transformation from the fixed to eigenframe is non-linear. The Gaussian nature is only 
preserved for u,xi,X2. In the Hessian eigenframe X12 — 0. From equations ( 16|17 1 



XiX2{xii — X22) = a:ia::2(Ai — A2) = 2x\X2W , 



(M + 2w)2i2 + (u-2wYx\ 



(24) 



In equation ( |15| l the averaging is now carried over the distribution of the eigenvalues with the measure 7r(Ai — A2) ( |Doroshkevich| 
1970) that accounts for eigenvalues being sorted, Ai ^ A2: 



8-2-7r 



dv (27r)Vl-7^ 



exp [-J7^/2] 



(Ai — \2)d\\d\2dx\dx2 Vs 5d(s) exp 



2(1-72) 



X2 



(25) 



^ here we made a choice of sign. Now in the coordinate frame that has the first direction aligned with the gradient of the field, i.e. 
a;2 = 0, to = (x\\ — a;22)/2, while in the frame that has gradient aligned with the second direction, x\ = 0, iS = (2:22 — a;'ii)/2 
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or in terms of u, w 

dC _ 



16 



exp [-77 /2I 



\w\dudwdx1dx2 Vs <5d(s)cxp 



2(1-7^) 



-2 

2^2 



(26) 



In the argument of the delta-function in equation ( 26 1 w can be zero only at special field points 
on a skeleton. So vanishing s requires either 



not at a generic point 
which describes, as expected, that in the Hessian eigenframe 



or X2 

one of the component of the gradient vanishes on the critical line. Since we have already chosen the coordinates so that the 



direction ' 
add up|^ 



1" is aligned with the largest eigenvalue and the critical lines can go in both eigen-directions, these two possibilities 



dC 
drj 



4^2 



(27r)3/Vl-7^ 



exp [—77^/2] / dww / dft (|2u; + {i| + |2i(} — exp 



2(1-72) 



- Aw' 



(27) 



Note that \w — u/2\ = |Ai| = |iii| and \u/2 + ™| = IA2I = |i22|. That is, the length of the critical lines per unit volume is 
given by the average absolute value of the Gaussian curvature of the field in the space orthogonal to the skeleton, given that 
in stiff approximation the direction of the skeleton is assumed to coincide with the gradient of the field. The reason for this 
is clear - the higher the curvature, the closer the next neighbouring segment of the skeleton can be, thus increasing the flux 
i.e. the length per unit volume. If we replace w —w in the second integral, we return to the formula (23 1 with integration 
over both positive and negative w. The integrated length of the critical lines is reduced to 



2V2 



wdw / du {\2w + u\ + 1211) — u\) cxp 



2 



(28) 



equations ( 27 1 and ( 28 1 are the results of the stiff approximation for the threshold dependent differential and the integrated 



lengths of the critical lines in 2D respectively. 



3.2 Primary critical lines in 2D: Skeleton and anti-Skeleton. 



The local skeleton is the subset of all the critical lines, which includes the parts that appear as the ridges in the field profile, 
rather than the valleys. This subset is described by the constraints that the skeleton lines should go along the largest eigenvalue 
Ai and, in addition, that this direction has the smallest curvature, |Ai| ^ |A2|. The anti-skeleton is a mirror structure describing 
the valley of the field and in all the results can be obtained by replacing rj —rj in the formulae for the skeleton. 

To derive the expression for the skeleton differential length let us return to equation ( |27[ ). The critical lines with Vp aligned 
with the largest eigenvalue direction have £2 = 0. Thus, only one term is selected by the fo-function: it is oc 2|A2| = 12™ + u\. 
The magnitude restrictions translates into u ^ 0, thus 



4^2 



dv (27r)a/2yr3:7 



:exp 



[—77^/2] / dww / du {2w + u) exp 



2(1-7^) . 



(29) 



This result should not be confused with equation ( 28 1, where w is integrated over full range of negative and positive values and 



which is strictly equivalent to equation ( 27 1, counting critical lines aligned both with the lowest and the largest eigen-directions. 



Performing the last two integrals one obtains for the differential length in closed form 



drj 



exp 



2 



1 + 



--I'll 



1 + Erf 



2sJ\-i' 



2x/27r 



■ exp 



2 2 

7 'n 



2(1 



and for the integrated skeleton lengtlj^ 



1 V2 

- + — 
8 47r 



0.23754 {y.R~^) 



(30) 



(31) 



Note that modulo the stiff approximation, equation (31 1 gives a universal, spectral parameter independent, scaling. Figure [s] 
demonstrates the threshold behaviour of the differential lengths for several values of the spectral parameter 7. 

The most important and robust result of our theory is the behaviour of the differential length at high density thresholds 



drj 



exp 



77? 



(32) 



It represents a bias similar to the one found in [Kaiser] ( |1984[ ) for the clustering of high critical points - maxima. According to 
the latter, the number density of peaks in regions above high thresholds is higher than on average. Similarly, the length density 
of critical lines above high threshold is enhanced relative to the mean. From the point of view of measurements, perhaps a 



^ If we do not sort the eigenvalues and, thus, do not restrict the w to be non-negative, then the notions of first and second direction are 
undefined, and we could choose now that the skeleton goes in, say, the first direction and X2 = 0. We will loose here factor of two which 
is recovered by having to extend w integration to negative values 

® In other words, one expect to find one segment of skeleton per linear section of (4.2iJ«). 
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Figure 3. left: dC^^^^ / drj / P{ri) (dashed) and g^skcl+antiskcl^gj^^pj-^-j (gQijj) \^ 2D for the spectral parameter values 7 = 0.3,0.6,0.95. 
Right: dC/drj/ P{ri) (solid) and its asymptotic behaviour (dashed) in 2D for the same spectral parameter values 



more interesting quantity than the differential length is the length per unit volume within the regions of high excursions of 
the field C{> rf). In terms of the cumulative length given by equation ([9|, C{> rj) — L — JZ{ri). Its asymptotic behaviour at 
high rj for the skeleton is found by direct integration of equation ( 32 1 



rErfc 



--1V 



(33) 



The first factor here is the fractional volume occupied by these high excursions of the field. Note that, at large rj the differential 
length divided by the PDF scales like rj'y/Rt = tj/Rq once the proper scaling with 1/_R, is introduced. Hence the differential 
length as a function of 77 together with the total length give access to two characteristic scales Rq and R^. See Appendix \K\ 
for a general proof of this result in N dimensions. 

The threshold dependence of the statistics of critical lines in the stiff approximation is determined solely by the spectral 
parameter 7. In the limit 7 = 0, when the distribution of the second derivatives of the field p is completely independent on the 
threshold, the length of the skeleton per unit volume within the regions with p/cro in the interval r],ri + dri is just proportional 
to the fraction of the unit volume that these regions occupy. Completely generally, for any type of critical line. 



dC 
drj 



(7 = 0) = 



exp 



2 



(34) 



When 7^1 the trace of the Hessian u becomes uniquely determined by the field level rj (recall equation (|6|). For over-dense 
regions with positive 77 equation ( |32[ ) is exact for 7 = 1, while no skeleton exists in under-dense regions in this limit. 
Near zero (mean density) threshold the dependence of dC^°^ /dr) is 



exp 



2 



1 VVT 

2 TT 



-72 1 
-— + 



-r 



-77? 



(35) 



Its details, in particular a step- like cutoff at negative r\ when 7^1, are sensitive to the definition of the primary lines. In 
under-dense regions with large negative densities the skeleton is exponentially suppressed. 

-77 for anti-skeleton, we obtain for the union of both primary critical lines 



Starting from equation ( 30 1 with 77 

^^skcl+antiskcl 



1 



with twice the integrated length 



exp 







1 


2 




4 



+ 



r skcl+antiskcl 



1 



= Erf 



777 + 



-r 



1 1 

4 V27r 



V2/r^ 

0.47508 (xi?~^) 



exp 



7^^ 
2(1-7^) 



(36) 



(37) 



This function is now symmetric in 77 with the skeleton providing the dominant contribution described by equation (321 in 
over-dense regions of space, and the anti-skeleton dominating the under-dense regions. Near the mean, zero, threshold of the 
field, both critical lines are present 



^^skel+antiskcl ^ -j^ 
977 



exp 



2 



2 2 
7 77 



1 \/r^ 

4 ^/2■K 2^2710"^^ 



(38) 



3.3 Secondary critical lines in 2D 

Secondary critical lines do not allow for a full analytical treatment and are investigated in Appendix [B] They are particularly 
important near zero threshold, since at this transitional regime the exact behaviour of primary or secondary lines depends 
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significantly on our somewhat arbitrary separation of tfie critical lines in types. In this paper we are tracking the skeleton 
— density ridges — as primary lines emanating from the maxima, until the largest eigenvalue ceases to be the shallowest. 
Alternative definition may, for example, somewhat extend the skeleton at the expense of secondary lines at lower densities 
as long as all the eigenvalues transverse to the gradient are negative, i.e until A2 becomes positive. As an advantage, the 
differential length of the skeleton and the corresponding secondary lines defined this way would not exhibit inflections at low 
densities that can be seen in Figures [S] and [BT] for high 7's. But the downside is that then one looses the ability to describe 
the primary lines analytically in a closed form. At the high density excursions the properties of the skeleton remain robust 
with respect to the variations in their exact definition. 

However the important advantage of the definition of the primary lines adopted in this paper lies deeper. The magnitude 
of the eigenvalue along the direction transverse to the gradient is connected to the stability of these trajectories near the 
critical lines and to their possible bifurcations. This is discussed in part in Section [5. 3[ 

Let us summarize the results for the total set of critical lines, primary and secondary combined, which are, of course, 
universal whatever the definition of the separate types. Summing up the results of this Section with the corresponding ones 
in Appendix [B] 



\/2 + acot(\/2) 



drj 

dC 
drj 



exp 



exp 



2 

!Z! 
2 



0.646071 {xR-^). 



V2(l-7^) + acot(V2(l-7^)) ^2(1-72) 



7r(3 - 272) 



(39) 
(40) 

(41) 



The full behaviour of the total differential length is presented in Figure [3] One should note the linear asymptotic behaviour 
at high density levels and the regular quadratic behaviour near zero density thresholcj^ Finally recall that in the section [s] we 
have omitted almost everywhere a 1/i?, factor for the quoted lengths and differential lengths. 



4 CRITICAL LINES OF 3D FIELDS 

In three dimensions, we carry the computations directly in the eigenframe of the Hessian, following closely the derivation of 
Sections |3.1| and |3.2| We present the formalism first for all the critical lines and then narrow our focus to the primary ones. 



4.1 The length of the critical lines of 3D fields 

In 3D, let us use the variables u — — (Ai + A2 + A3), to = (Ai — A3)/2,t; — (2A2 — Ai — A3)/2. In the Hessian eigenframe 

= (A2 — A3)i2i3 = (to + v)S:2X3 , = (A3 — Ai)a;ii3 — —2wx\xz , — (Ai — A2)xi£2 = (to — v)x\X2 , (43) 



and 

— 1 
Vs 

— 2 
Vs 



{0, A2(A2-A3)i3, A3(A2 - A3)£2} = <{ 0, ~ ^{u - 2v){'w + v)xs, ~ ^{u + v + 3w){w + v)x2 



{Ai(A3 - Ai)i3, 0, A3(A3-Ai)ii} 
{Ai(Ai - A2)i2, A2(Ai-A2)ii, 0} 



2.2 

-{u + V — 3w)wx3, 0, -{u + V + 3w)wxi 



^{u + V — 3to)(to — v)x2, 



^{u — 2v){w — v)x3, 



(44) 



In the eigenvalue space the measure is 27r^|(Ai — A2)(A2 — A3)(A3 — Ai)| and the eigenvalues are considered sorted. For sorted 
eigenvalues the choice of the directions has been fixed and the (s'^, s^), (s^, s"^) and (s^, s^) pairs of surfaces describe different 
possibilities for the critical line. Those choices add together in the average integrated length. Using the variable TO,i3 the 
condition of eigenvalues being sorted is to ^ 0, —w ^ v ^ w. 

Let us consider the critical lines that are the intersections of (s^, s"^). Their differential length is given by 

dC „ 2 3 S^/^IS^S'/^ 

-— — 2iY ■ - ■ , exp 

dri 2 (27r)Vl-7^ 



1 2 
o»7 



(Ai — A2)(A2 — A3)(A3 — Ai)| dAidA2dA3da;idi2da;3 Vs x Vs (5D(s^)fo(s^) 



X exp 



(u - -yri)^ 15 .2 5 .2 



3-2 3 _2 3 _2 

-a;, Xy x-i 

2 2 2 



(45) 



For all 7 but 7 = 1, for which 

1^(7 - l,r, ^ 0) ~ ^ exp [-„V2] (I + :r^\r,\' - -^\vf ■ ■ ) 
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Figure 4. Left; The skeleton dC^^°^ /d'q/ P(ri) (solid) and its asymptotic behaviour at high density thresholds (dashed) in 3D. The 
anti-skeleton is described by the curves symmetric with respect to a reflection of r). Right: dC/drj/ P{ri) its asymptotic behaviour for the 
total set of critical lines. The spectral parameter values are (from bottom to top at high rj) •y = 0.3, 0.6, 0.95. 



Integration over (Sd(s2) and 60(33) leads to the only possibility £2 = 0,xs — 0. That is, the choice of the surface 5*2 and S3 in 
the Hessian eigenframe describes the skeleton along which the gradient is aligned with the direction 1, correspondent to the 
largest eigenvalue, while in the directions 2 and 3 the components of the gradient of the field vanish. With X2 = X3 = we 
get a simple expression for 



j Vs2 X Vssj = 1A2A3(A3 — Ai)(Ai — A2)| xl = -\{u — 2v)(u + v + 3w){w ~ v)w\ x\ , 
while the subsequent integration over X2 and xz using fo-functions and afterwards over x\ gives 



7r(l — 7^) 



cxp 



1 2 



dAi dA2dA3 ( Ai - A2 ) ( A2 - A3 ) ( A3 - Ai ) 1 A2 A3 1 exp 



(u - "jvif _ 15 .2 
'2(1-7^) 2™ ■ 



5-2 

-V 

2 



(46) 



(47) 



Notice again that what the integrand involves the Gaussian curvature in the direction orthogonal to the gradient, which in 
stiff approximation is the direction of the filament itself. The contributions of the critical lines directed along the second and 
third eigen-direction is given by similar considerations and are added together when all critical lines are considered. Changing 
variables one finally obtains 



dC _ 3^5^/2 exp [77V2] 



dr) 47rV27r(l - 72) 



/ 



dudwdv w{'w — u ) (IA2A3I + IA1A3I + IA1A2I) exp 



(a - 777)^ _ 15 -2 
2(1-7^) 2™ 



5-2 

-V 

2 



while the integrated length is 
3355/2 



L = 



47r2 



j dudwdv wiyo^ — v'^) (IA2A3I + IA1A3I + | A1A2I) cxp 



1-2 15-2 

u w 

2 2 



5-2 

-V 

2 



withd] 



IA2A3I + IA1A3I + IA1A2I = - [|(u - 2i}){u + i + 3w)\ + \{u + v + 3w){u + i- 3i5)l + \{u - 2{))(fi + S - 3w)|] 



(48) 



(49) 



(50) 



The equations (481 and (49 1 account for all the critical lines. In Figure |4] (right panel) the results for 3D critical lines are 
plotted while the discussion of the corresponding asymptotics is given in the Appendix [C] We shall now turn our attention 
to the study of the primary lines and, in particular, the 3D skeleton that delineates the over dense filamentary structure and 
is of more direct observational interest. 



4.2 Primary critical lines of 3D fields: Skeleton and Anti-Skeleton 

The subset of critical lines identified with the skeleton correspond to the lines with the gradient aligned with the largest 



eigenvalue Ai while having Ai + A2 ^ 0. In equation (481 such lines are described by the first term ~ jA2A3j. The differential 



One should note the correspondence with the well-known result for the number density of extrema of the field i jBardeen et al.|l986| 

5 . 



-/Voxt oc J dudwdv w{w'^ — v^)\XiX2Xz\ exp 
which is determined by the mean three-dimensional Gaussian curvature IA1A2A3I 



1-2 15 _2 

u w 

2 2 
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w 



1.0 






\ ^^^^ 






\ 0.6 






\ 0.4 






\ 0.2 







-0.6 -0.4 -0.2 0.0 0.2 0.4 



Figure 5. Left: Integration zones in the v, u plane for the 3D skeleton analysis. Variables arc given in units of w. v varies from —ui to 
+w, while u must be greater than ^{v + 3w). In the allowed shaded region > A2 ^ A3 everywhere. Horizontal lines mark the further 
subdivision of the integration space if the order of integration is changed according to equation | |53| l. Right: Integration zones in the 
{v, w) plane after ii has been mapped to the [0 — 00] interval. Variables are given in the units of u. The lower triangular zone corresponds 
to the semi-open rectangular band above the red dashed line in the left panel. In this region the integrand is given by the first term of 
equation \53\. It dominates the high rj asymptotics. 



length of the skeleton is then 

g^skel ^ 3355/2 ^_^2^2] 

dri ^ 47r2 ^27r(l - 7^) 



3-5''^^exp [-7772] 



47rV27i-(l - 7^) Jo 



roo rw poo 

I dw dv du 

Jo J -w J ^{v + 3w) 

/W POO 
dv / 



wiw^ — t;'^)A2A3 exp 



(m - in? 15-2 5-2 



du with —V ) (m — 2«) (m + -C + 3t5) exp 

i{5+3iB) 



(u - jri) 15 .2 5 .2 

-2(w) -y" --2\ 



(51) 



The integration in v~u plane is limited to the region u > |(w + 3w), as shown in the left panel of Figure [s] The integrated 
length of the skeleton is 



T^skei ^ 0.046186 {xR~^) 



(52) 

that is, one expect on average one skeleton line crossing a random ~ {5Rt? surface element. The results of integration of 
equation (51 1 are presented in the left panel of Figure [4] 



4-2.1 Asymptotic behaviour at 777 ^ 00 

To study high rj asymptotes it is useful to change the order of integration to have the u integral as the outmost one. The 
inner integration in v-w plane is then carried out over the region shown in the right panel of Figure [5] 



roo ru! roc roo /'^/2 ru! roo f-u r2u — 'd'W 

I dw I dv I du ~* I du I dw I dv + I du I dw I dv 

Jo J -w J h(v+'iu>) Jo Jo J -w Jo Ju/2 J- 



10 J -w J ^(v + 'iw) 

The last term is exponentially suppressed as 77 — + cxj while the first one gives 
dC"""' j^-.oo 3 • 55/' exp [-r;V2] 



(53) 



47rV27i-(l - 72) 



roo roc rw 

du dw / dv w{uP — iP'^iu — 2v){u + u + 'iw) 

Jo Jo J -vj 



exp 



{u - 7?7) 15 ,2 5 .2 



exp 



1 2 



(yvf + 9i'yri)/VT0^ + (9/10 - 7') 



(54) 



The leading quadratic and the next linear terms can be recovered found by replacing u — > 777 in the pre-exponential factor 
and treating the exponent as the ^o-function. A more detailed asymptotic study of this Laplace-type integral is required to 
recover the third-order constant term, that also contributes to the accuracy of the expansion at the level demonstrated in 
Figure |4] 

One finds that in the leading order in 77 the skeleton has the differential length growing as (777)^ (see also Appendix[A| and 
involves, as expected, a third of all the critical lines (compare with Appendix[c| in the regions of high excursions concentrated 
around the maxima of the field. However, at intermediated thresholds, the skeleton constitutes more than a half of all critical 
lines, highlighting enhanced importance of the filamentary dense ridges among other critical lines. 



Note the appearance of the linear in 77; term in the next to leading order for the skeleton, that canceled out for the critical lines. 
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Figure 6. First coefficients of low-r; power expansion {Aq - blue, A2 - yellow, A4 — green) of the differential lengths of the 3D skeleton 
{left), inter-skeleton (middle) and total set of primary critical lines (right). The odd terms (e.g. Ai — dashed) are present only for the 
asymmetric case of the primary lines. 



4-2.2 Power series at r] ^ Q and Hermite expansion 



Using two alternative series representations of the shifted Gaussian form that encodes the dependence of the skeleton on the 
threshold rj 



exp 



(?i - 777)^ 



2(1-7^) 



27r 



exp 



u 
Y 



: exp 



2(1-7^) 



°° 1 
5Z ~7f 



V2^(l-7^) 

we obtain either power series or HermitfQ ( [Novikov et al.|2006[ ) expansion of the differential length 

o/-skGl 1 

- = ^exp [V/2] X 



drj 



(55) 
(56) 

(57) 



where 
Aki-y) 

and 



3 • 5^/2 



4-7^2^/fc!(l — 7^) ''2^ Jo 



dvj / dv 



du w('w — V )(u — 2v)(u + V ~\- ■iw) 



X exp 



2(1 - 72 



15 _2 



5-2 



3-5 



5/2 



dw / dv 



du 



du w{uP 



V ) (m — 2«) (« + « + 'iw) exp 



u 15 _2 5 .2 

W V 

2 2 2 



Hk{u) 



(58) 



(59) 



These two expansions are similar but distinct. The power-law expansion is suitable for an accurate analysis of the differential 
length near zero threshold for all 7 < 1. On the other hand, the expansion in orthogonal Hermite polynomials is useful as an 
approximation over an extended range of thresholds. Both series are improper for 7 = 1. 

Although these coefficients can be computed analytically, their expressions are too cumbersome. Instead, we plot several 
leading ones in Figure[6] Remarkably, the power in Hermite expansion is concentrated in a few low order terms, in particular, 
k = 0,1,2,3 for the skeleton, with subsequ ent terms forming a slowly decaying oscillating series. This finding confirms in 
3D the conjecture of 



Novikov et al. 



(20061. The contribution of the first three most dominant terms, "^^^q BkHk{ri) 



0.0462 -|- 0.07 51777 -I- 0.04647^ (r;^ — 1) has the same structure and remarkably similar coefficients as the high rj asymptotics of 
equation (54 1 which evaluates to 0.0477 -I- 0.085277/ -I- 0.05317^ (t;'^ — 1). This explains why the high 77 asymptotics provides a 
visually good fit through all thresholds when 7 is not too high. At 7 ^ 1, the oscillatory tail of Hermite series provides the 
correction that refiects the irregular nature of the expansion in this limit. 

The power series expansion refiects the features of the Hermite expansion. Starting, by definition, at ^o(O) — Bo — 
jCl^t = 0.0462, Aq beh aves as ^0(7) « rtot(l - 7^) over most of the 7 range. Coupled with ^2(7) « const = jCtot and 
Ai(7) ^ 0.0751^1-7^ we get for the first three orders ELo Akiy){-fvf ~ 0.0462 + 0.0751^1 - 7^7?? + 0.04627^(772 - 1), 
close both to the Hermite expansion and to the high 77 law for moderate 7. On the other hand, the power series expansion 
explicitly demonstrates the increasing importance of higher-order terms for 7 > 0.8. 



^ we use here the normalized Hermite polynomials following probabilistic definition, l/\/27r ducxp [— l"'^] Hk(u)Hm(u) = 5^^. 
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Figure 7. Differential length dC/d'q/P(r^) of the intermediate (Jeft panel) and combined primary lines {right panel) as function of the 
threshold rj in 3D. Different curves from blue to yellow correspond to the spectral parameter values 7 = 0.3, 0.6, 0.95. The dashed curves, 
drawn only for positive rj, correspond to high-r; asymptotic solutions. 



4.3 Primary critical lines of 3D fields: Inter- Skeleton and the overall behaviour 

The intermediate primary critical lines are associated with saddle-like regions where the largest eigenvalues in magnitude are 
Ai ^ and A3 0, and have opposite signs, and the shallowest direction aligned with the gradient is the second one with 
— Ai < A2 < —A3. Their appearance reflects the complexity of critical lines in space of more than two dimensions. 

The differential length of the intra-skeleton computed in the stiff approximation is presented in Figure [7] The conditions 
for intermediate lines are prevalent for the regions of the field of moderate values - within 2a {\ri\ < 2 of the zero mean for 
7 > 0.6. Although the occurrence of the intra-skeleton within these regions is never large {dC/dr]/ P{rj) is relatively small), 
the regions corresponding to a near mean density occupy large fractions of the total volume, and as the result the total length 
of the intermediate skeleton is almost twice that of the skeleton or the anti-skeleton: 

jr^intcr ^ 0.087533 {xR''^) . (60) 
It constitutes nearly a half of the total length of the primary critical lines 

^prim ^ ^skel ^ ^antiskol ^ ^intcr ^ 0.179905 {X R^'^) . (61) 

At high \r]\ thresholds, in very dense regions near maxima or under-dense regions near minima of the field, the intermediate 
skeleton is rare. 

The total set of the primary critical line is even more than the skeleton dominated by the low order terms in Hermite 
expansion. Indeed, Figure |6] demonstrates that just the first two terms (odd orders are absent due to symmetry) in Hermite 
series are dominant, YlkLo BkHkiv) ~ ^tot^ (l + 0.3407^(77^ — 1)). 



4.4 Validity of the stiff approximation 

Let us consider the opposite to "stiff" regime, when the derivatives of the Hessian dominate the Vs, 



(VmS' 



7 



-1 \ ^ ij _ 



(62) 



jki 



Although not natural for cosmology-inspired spectra, such a situation arises when the power spectrum has an extended short 
wave tail with spectral inde^^ n between —9 and —5. Such spectra have small 7, _R ^ i?, and there are many infiection 
points of the field per extremum. Interestingly, this regime also automatically means that the correlation between the gradient 
and third derivatives of the field is small. 

Using the Hessian eigenframe formalism, we can obtain the important results without explicit computation of the differ- 



ential length. Let us focus on the critical lines corresponding to the first eigenvalue. Equations (431 for 5-surfaces gives rise 
to two fo-functions, Su{'2wxiX3)6n{{w — v)S:iX2) that after integration over the transverse gradient components £2 and X3 
enforce X2 — X3 — 0, with the Jacobian factor l/\2w{w — v)xl\. The length element in this frame obeys 



E 



2il 3jl 



: 7 '^xtll>{Xklm) ; 



(63) 



In a cosmological framework this takes place when the density field with n < — 1 spectrum is smoothed with a top-hat window. 
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where the last expression defines the ipi^Xkim.) function. The differential length is now given by 



1 



3355/2 



: exp 



1 2 



dudwdv (w + v) exp 



[u - 777)^ 15 -2 5 .2 

, TT^ W V 

2(1-72) 2 2 



X-i^dXld^ Xklni1p{xklm)Pl{xi,X2 = Xz =Q,Xklm) 



(64) 



The last integral, with Pi given by equation (D6l, is a function of 7 only. The first term shows that, since the integrand 
prefactor is independent on u, the differential length does not depend on the threshold 77 at large 77 (it does at small r) only 
because of non-trivial integration boundaries dependent on the exact type of critical lines). This is not surprising, since in 
this limit, there is little link between the skeleton length and the second derivatives, the only ones that are correlated with 
the field value. Such threshold independent behaviour is not observed in simulations with cosmological spectra, which argues 
once again for the statistical validity of the "stifT" approximation. 



4.5 Measurements 

In this section we compare the predictions of the local theory in stiff approximation with the measurements of the statistical 
properties of the critical lines done on realizations of the Gaussian fields with different power spectra. 

We perform the measurements on critical lines found according to the global definition. The measurements are carried as 
follows: a set (typically ~ 100) of scale-invariant Gaussian random field of a N-D maps (typically 10242) or cubes (typically 
256^) is generated with a given power index of n = 0, —1 or —2. The N-D cube is then smoothed via convolution with a 
Gaussian kernel of width 6 pixels. The spectral parameters, 7, 7 etc... are computed through the second moments of the 
derivative of the smoothed field. The set of critical lines is then extracted as the intersection of the peak patches and void 
patches (see Sousbie et al. (20081 for details). In Figure [s] an example realization of the primary critical lines in 3D cube is 



shown. Since the algorithm produces a set of segments describing those critical lines tagged by the underlying (smoothed) 
density field, it is straightforward to compute the total and differential length per unit volume of the whole set. The differential 
length per unit modulus gradient is extracted by tagging the critical lines with this modulus (obtained via Fourier transform 
differentiation) and proceeding as before. Finally, the curvature of the skeleton is measured by computing the local curvature 
of a set of adjacent segments via finite difference. 

Let us emphasize that these measurements correspond to properties of the global skeleton, whereas the theory developed in 
this paper is focused on the local skeleton. Hence even more remarkable is the match between the measured and the theoretical 
differential lengths for all values of 7, that is exhibited in Figure [9] This accuracy should be considered as indicative of the 
correspondence between the stiff approximation to the local theory and the global set of critical lines. 



5 OTHER STATISTICS AND SPECTRAL PARAMETERS 

In the previous sections, the emphasis has been on the differential length of the critical lines as a function of the excursion in 
density. As argued in Sousbie et al. (20081 and demonstrated here, it provides means of constraining the shape parameter, 7. 
Let us now explore other statistics which will allow us to constraint other shape parameters. In particular, let us demonstrate 
that the differential length as a function of the excursion in the modulus of the gradient of the density and the differential 
curvature depend on the second shape parameter, 7. Finally, we investigate the number density of singular points on the 
critical lines. 



5.1 Differential length versus the gradient modulus 



The differential length of the skeleton with respect to the threshold 77 carries information on the spectral parameter 7 thanks 
to the correlation between the value rj and the Hessian of the field. In the stiff approximation the Hessian curvature completely 
determines the length of the critical lines. For the exact formulation, the length also depends on the third derivatives, that are 
correlated with the first derivatives via the parameter 7. Thus, measuring length as a function of the modulus of the gradient 
should carry information on 7 and provide an estimate of an impact the third derivatives have on the length statistics of the 
critical lines. 

To demonstrate the dependence of the skeleton length on the gradient of the field in "stiff" approximation let us return to 



equation (45 1 which we take integrated over all density thresholds. As before, we perform the integration over the 5-functions 
that enforces alignment of the gradient with the first eigen-direction, X2 = x^ = Q 
but rather take the differential of the result with respect to xi. Noting that = 



however this time we do not integrate over 

: X 



equation (47 1 



x'^ + X2 + Xg. we obtain in place of 



dC 3 3^''^15^5^''^ 



exp 



2 



(Al - A2)(A2 - A3)(A3 - Ai)| dAidA2dA3|A2A3| exp 



1-2 15-2 5-2 

It W V 

2 2 2 



(65) 



dX 2 (27r)5/2 

where the last integral does not depend on X. Dividing by the integrated length, L — dC/ dXdX , and generalizing the 
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Figure 8. An example of set of primary critical lines (resp. skeleton in blue, intermediate in magenta and anti skeleton in gold) for a 
scale invariant power spectrum with 7 = 0.6 in a 256^ box smoothed over 5 pixels. 



dL/d77/P(77) dL/d77/P(77) 




Figure 9. The relative differential length, dC/dri/PDF, measured in simulation of 2D (Jeft) and 3D (right) Gaussian random fields with 
scale invariant power-law spectra versus predictions of the local theory in stiff approximation (solid curves). The spectral parameter 
7 = 0.71, 0.59, 0.39 for the 2D and 7 = 0.77, 0.70, 0.60 for the 3D simulations. 



result to fields in arbitrary A'^ dimensions we conclude tliat 

stiff 



Ldx) 



cxp 



2"^ 



(66) 



The exact dependence of the differential lengtiis will deviate from tliis form in a 7 -dependent way. It is natural to 
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Figure 10. left: measured dC/dX/L as a function of X = | Vp| for 7 = 0.71, 0.58 and 0.38 using the set of 25 2D simulations of 1024-^ 
Gaussian random fields with scale invariant spectra smoothed over 7 pixels. Right: value of the fit parameters (see equation (j67|). 
Note that cq = 1. 



parameterize such deviation expanding the true statistics in Hermite series around the stiflF approximation 



1 dC 
LdX 



exp 



C2fc(7) 



(67) 



This choice of expansion is dictated by the orthogonahty of the Hermite polynomials with the weight oc exp[— \/iVX^/2] on 
the interval X ^ 0. Thus, cq = 1. If the deviation from the stiff approximation is small, one expects the expansion to be 
dominated by the n — term, while the subsequent terms should quickly fall in a orderly fashion. 

To gain understanding on how the coefficients C2fe(7) behave with 7, let us consider again the lax situation, opposite to 
the stiff case, when the third derivatives of the field dominate the length statistics. Our starting point is equation ( 64 1 which 
has the following structure when we consider the differential length with respect to the = X 



— OCX exp 



3 -^2 

2^ 



dui < 



exp 



3 (^^1 - 7X)= 



+ exp 



3 {u^+jXf ^ 
2 1-72 



k1^'{Xijk)Pl{x2 = X3 = 0,Xijk/ui) 



(68) 

where Pi is given by equation (D6I with the dependence on ui factored out. The difference with the stiff approximation is 
large even for 7 = as the gradient's dependence becomes oc exp[— 3/2X^] in place of the stiff scaling oc exp[— 3/2X^]. 
Using now this factor as the weight, for 7 7^ we expand the expression in the brackets in generalized Laguerre polynomials. 
The expansion coefficients are of the form 7'^* exp[— 3ui/2] X]m=o dm^^"^ H2m{V3ui); denoting the result of the integration of 
the expansion coefficients and all of the residual factors over the third derivatives by *I'fe(7) we obtain 



LdX =3V^^^"P 



3 -^2 
2^ 



E 



(2fc + 1)!! 



(69) 



where, again, ^'0(7) ~ 1- With the help of the relation between the Laguerre and Hermite polynomials 
3X\\Li^^^^ (3/2X') = (-l)''-2-'= (//2fc+2(\^X) + (2fc + l)i/2fc(V3X)) , 



we can cast equation ( 69 1 in the form of equation ( 67 \ 

I 

(2fc + l) 

3 



1 ar'"" 

LdX 



exp 



3 ^2 
'2^ 



E .ai + i)!! ^'' {H2k+2{V3X) + {2k + l)H2k{VsX)) *,(7) , 



exp 



-X' 



1 + E L^-lTi! [7''"'**-i(7) - 7''*fc(7)] H2k{V3X) 



^ (2k- 
fe=i ^ 

The coefficients 02^(7) are 

CO = 1 , C2 = V2 (1 - 7'*i(7)) , C2k = ((-l)'=-'/(2fc)!/(2fc - 1)!!) f*'^ (*ft-i(7) - 7'*fe(7)) ■ 

In particular, in the limit 7-^0 the first two coefficients remain finite and of equally significant magnitude co = 1, C2 
while all the other ones vanish. 



(70) 

(71) 

--V2, 



Figure 10 and Figure 11 present the measurements of dLjdX jL in 2 and 3D respectively, together with the corresponding 



coefficients given by equation ( |67[ ). It is found that these coefficients are significantly smaller in two dimensions, a clear 
indication that the stiff approximation holds better in 2D. 
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Figure 11. left: measured dC/dX/ L as a function oi X = |Vp| for 7 = 0.86,0.83 and 0.79 using the set of 3D simulation of 128^ 
Gaussian random fields with scale invariant spectra smoothed over 7 pixels. Right: value of the fit parameters c^., Note its faster 
convergence combined with a larger amplitude relative to the 2D case. 



5.2 Statistics of the curvature of the critical lines 

The local curvature, k, at a point on a curve specified by the tangent vector u — dr/dt is determined by the acceleration of 
the tangent vector u = du/di = u • (du/dr) transverse to the curve direction: 

^ ^ |u X u| ^ |u X ((Vu) ■ u) I 
up' |u|^ ' 

Importantly, the curvature does not depend on parameterization t, nor on normalization of the tangent vector u. In the local 
theory, the tangent vector to a critical line is orthogonal to Vs'(a;fc, Xki) and can be taken to be 

u = e ■ Vs (2D) , u = Vs^ • e ■ Vs^' = Vs' x Vs-' (3D) ; (73) 

so the curvature k is the random quantity which involves the derivatives of the field up to fourth order, 

(2D) K = 1V.S- (VVs)- Vsl/IVif , (74) 

(3D) «: = |(^Vs* X Vs^) X [(^Vs* X Vs^) • V (^Vs* X Vs^)]|/|v.s* X V.s^p . (75) 

The curvature of the critical lines fundamentally refiects the derivatives of the field higher than the second. If they are 
neglected, the curvature is identically zero. Explicitly, the contributions that do not involve higher derivatives, in 2D 

(2D) (uxii)'' = 4xlxlxt{Xi~\2fxt + ..., (76) 
(3D) (uxii)' = {xlxl + xlxl) {\i - X2)^ (Xi - Xs)^ Xl {3x1x1x1 + {xlxl + xlxl) XiY + ... , (77) 

vanish when the correspondent critical line conditions a;2 = or 3:2 = 2:3 = are appliecj^ 

The integrated curvature over the length of the line, C — J ndL is a useful dimensionless characteristics of the overall 
extend a line is curved. We have seen that the critical line length in volume dV is dL oc \\/s\S-D{s)dV and dL oc |Vs* x 
\/ S-' \Su{s^)S-d{s^ )dV in 2D and 3D respectively. Averaging over statistical distribution in regions above threshold r] we obtain 
the mean density of the integrated critical line curvature C = (dC/dV) 

(2D) C(J7>) = — / dxd'^Xkd'^Xkld'^Xklnid^XklmnK.{xk,Xkl,- ■ ■)\^s\3-D(s)P{x,Xk,Xki,- ■ ■) , (78) 



r]>x 



(3D) C(r7>) = 1^ dxd^Xkd'^Xkid^°Xkimd^''xkirnni^{xk,Xki,---)\\7s^x\/s^\s-Dis^)5-Dis^)P{x,Xk,Xki---), (79) 

where the integration is carried over all the derivatives up to the fourth order. The required joint probability function is given 
in equations (D19I and (D21I. 

Let us consider 2D case and estimate the curvature by using stiff approximation for the tangent vector it while following 
its local variation which involve higher derivatives of the underlying field. In the Hessian eigenframe, assuming the skeleton 
lies along 1, if u is approximated by its stiff counterpart we have: 

stiff 

= l(Ai + A2)xii2| , (80) 



Note that by construction, the torsion: r = |u • (u X ii)|/|u X lip contains only the terms proportional to at least the third derivatives 
of the field. 
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Figure 12. !/(/?*(«;)), the mean curvature radius in units of i?* as a function of r] measured in simulation of 2D (Jeft) and 3D {right) 
Gaussian random fields with scale invariant power-law spectra. The top curves correspond to spectra with more power at small scales 
and higher 7 and 7 spectral parameters. 



so that (taking into account the measure in the eigenframe and the So function of 5* in 2:2) 
dri 

, ^ / . „ „ . r H \ ^ 



TT / d\ld\2d Xijk\{Xi + X2)xii2\ Poiv,^kl)PliXijk) 



4ttR,R 



v/2(l-7^)cxp 



2 2 

7 7] 



2(1-72) 



777 



^2(1- 7^) 



777 



exp 



2 



(81) 



the last evaluation being done for the primary critical lines. 

We have measured the mean curvature of the skeleton lines at the threshold rj, 

_ dC/dr) 

= wm ' ^ ^ 

in simulations of the Gaussian random fields of different spectra, using the global skeleton techniques. Figure \12\ displays 
the 2D (left panel) and 3D (right panel) results in terms of the curvature radius, 7?curv = 1/{k). The measurements show 
that for the spectra we consider, the averaged curvature is insensitive to the density threshold for low-to-moderate threshold 
values showing a plateau in the interval — 2 < < 2. This indicates that in this regime the curvature of the skeleton does not 
depend on 7, but rather on 7 and perhaps 7. It follows that in 3D the critical lines are relatively more wiggly than in 2D. 
If we use the lower value of RfSr v a s a guidance, it seems the stiff approximation is less accurate in 3D that in 2D, as could 



be expect ed. The stiff estimate (81 1 gives the threshold-averaged mean density of the integrated curvature C 

as 



and, 



equation (37 1, the mean curvature radius, i?^ 



2nR,R 



He 



L 

C 



V2 + 7v/2_^ 



2.985- 



--R, 



using 



(83) 



V1-7V2 ■ V1-7V2 

This result captures the qualitative dependence on 7 observed in simulations, but is a factor of three smaller in the magnitude 
of the curvature radius. This shows that the global skeleton, used in the numerical measurement, is notably straighter than 
the local critical lines, although the dependence of curvature on the spectral parameters seems similar. 

Note in closing that in 2D, (resp. 3D) the knowledge of the differential length, curvature (resp. length, curvature and 
torsion) corresponds to an exhaustive global statistical description of the critical lines. 



5.3 Singular points of critical lines 

Let us now ask ourselves the following question: are there any special points along the skeleton? The obvious ones are the 
extrema of the field itself where critical lines intersect. Beyond this, one can anticipate two other types of singular points. 
The first type corresponds to points where the curvature transverse to the direction of the critical line vanishes along at 
least one axis: typically, in 2D, they mark regions where a crest becomes a trough, or vanish into a plateau. The second type 
correspond to points where the critical lines would split, even though the field does not go through an extremum: a bifurcation 
of the lines occurs along the slope; the occasional skier or mountaineer will be familiar with a crest line splitting in two, even 
though the gradient of the field has not vanished. From the point of view of the theory of random fields, the frequency of 
such points is an interesting venue: indeed we expec t th at steep power spectrum present relatively more bifurcation points 



as R, the distance between infiection points (see Sec. 2.21, becomes much shorter than Rt, the distance between extrema. In 
an astrophysical context, the statistical properties of the first type of points, and in particular their clustering properties are 
of interest for understanding the geometry of galactic infall, which in turn is believed to play an important role in defining 
the morphological properties of galaxies. The multiplicity of the maxima (i.e. the number of connected skeleton segments) is 
also of interest in the context of galaxy formation and feedback. In more abstract spaces, such as position-time, identifying 



bifurcations is important to pin down merging events (see e.g. Hanami (20011 and Appendix [A| 
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5.3.1 Defining the skeleton singular points 

Formally a singular condition along the skeleton occurs when at some point the determination of the critical line direction 
fails. It means that at this point the matrix Vfc5" of equation ^ has more than one distinct right null-vector, or, equivalently. 



all defined by equation ( A7l are zero. The only case when it happens exactly is at the extremal points of the field Vp = 0. 
There are no other formal singularities on the local critical lines, since when Vp 7^ 0, the requirement M*^ = sets A'' relations 
between the field gradient, second and third derivatives which have vanishing probability to be simultaneously satisfied along 
a line in a random field. 

The failure of the formal definition to identify all the physically interesting situations primarily refiects the inadequacy 
of the local skeleton construction, which only utilizes locally quadratic approximation to the field, to map the field near the 
singular points. [^Figure 13 gives a 2D example. In 2D, the critical lines are zero levels of the scalar S-function, while VS — 
at the extrema of S field. In Figure [13] the region where a ridge splits into two is shown. One expect two critical lines cross 
there, with three branches following the ridges, and one following the through between two of the split branches. Instead, the 
locally defined critical lines are not allowed two join at the bifurcation point since the formal condition VS* = is satisfied 
just off S — contour, rather they artificially reconnect near the bifurcation point into two non-intersecting segments. 

We conjecture that the critical lines experience a qualitative change in behaviour in the vicinity of the points where either 
the Hessian eigenvalue of the orthogonal to the gradient direction vanish, or becomes equal to the one along the gradient. 
Namely, if, for definiteness, Vp is taken to be along the first eigen-direction, A2 — 0, or A2 = Ai. We call the first case the 
"sloping plateau" as it designates the entering of a fiat region, and the second, tentatively, the "bifurcation" as it designates 
the places of possible reconnection of critical lines. In particular, at the A2 = Ai points most of the transitions from primary 
to secondary behaviour take place. Remarkably, these special points on the critical lines are recov ered by the formal singular 



condition \M | = if Vfc5" is evaluated in the stiff approximation. As given in equation (A12l, along the ND critical line 
defined by a;2 = . • . = = 0, IM'^''*^! = a;^~^ ni>i i^i ~ ^i) — gives rise to three classes of situations: (i) xi — 
corresponding to extremal points; (ii) one of Ai = corresponding to slopping flattened tubes; and (iii) one of Ai — Ai, 
corresponding to an isotropic bifurcation. 

Since it is beyond the scope of this paper to develop the full theory of these special points, we will focus here on their 
number density for isotropic 2D Gaussian random fields, leaving more detailed investigation to future work. 



5.3.2 Number density of the singular points of the 2D critical lines 

In 2D, the skeleton's singular points correspond to points where Sk = VfcS' = 0. The number density, nB(??) of singular points 
below the threshold rj is equal to 

"-b(?;) = / dxd'^Xkd^Xkid'^Xkimd'^XkimnP{x,Xk,Xki,- ■ ■)\<iet{'Vk'Vis)\5o{si)5o{s2) . (84) 

r]'>x 

The simplest case of the skeleton singular points VS = are, according to equation the extrema of the field itself, 
xi = X2 = 0. Indeed when both xi and X2 vanish 

|det(VfcV;s) |5d(si)5d(s2) = \xki\5i,{xi)5D{x2) , (85) 

which is exactly the integrand involved in the number density of extrema of the field. The extrema number densities, for 
reference, are given in 2D by ( |Longuet-Higgins|p57| 



: exp 



2^Vl - 2773 V 2(1-272/3) 



(86) 



9n-min+max Sngaddlc , 1 2^ , , 2x 1 I t a'7\ 



The singularity of the extrema from the points of view of the critical line theory is manifest in the fact that at extrema several 
critical lines intersect. 

The gradient of S, evaluated in the stiff approximation, in the Hessian eigenframe has the components 

sf'** — X2\i (Ai — A2) , and sf'^ — xiM (Ai — A2) . (88) 

and involves only second derivatives of the fleld. Remarkably, within this approximation, there are new singular points that 
lie on the (local) critical lines. The reason is that among two conditions needed for VS to vanish, one is already automatically 
satisfied by being on a critical line. 

To be specific, let us consider the critical line that corresponds to the X2 =G condition in the Hessian eigenframe. Then 
sf*'*^ vanishes everywhere along this line. The requirement sl''^ ~ has a solution at the extremal points, xi — 0, but also in 
two other cases, namely A2 = or A2 = Ai, that we conjectured to be of interest. 



^2 Similarly, the bifurcation points for the global fully-connected skeleton ( jSousbie et al.|2008^ also formally merge with critical points 
in tile strict sense due to sliarp topological theorems I Jost|2008 1. However they appear if the skeleton is viewed with a finite resolution. 
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Figure 13. Left: the three types of singular point on the critical lines (in solid blue: primary; dashed: secondary; green: gradient lines of 
the global skeleton): an extremum (open circle), a "bifurcation" (white square), and "slopping plateaux" {black squares). The thin lines 
are the isocontours of the field. Right: Detailed view of the bifurcation region: The pink and purple lines mark the conditions xn = X22 
and X12 = 0, which when intersect give the point where Ai = A2. This point is a singular point of the critical line (white square). The 
gold line is the condition A2 = that marks the "sloping plateau" on the critical lines. The red and blue dashed lines are zero isocontours 
of two components of VS"^''** . The VS"'*'''^ = criterium pin-points exactly all types of the singular points on the critical lines. The black 
cross marks the position of the point \7S = 0. 



The first situation, the sloping plateau with a flat transverse gradient, only occurs on secondary critical lines since it 
implies Ai > 0, and corresponds to 



det 



a::ia::ii2a;222 



(Ai - A2) 



(89) 



hence 



dr] 



-i- J dAlPo(77, Ai, A2 = 0) X j dXid^Xklm,Pl{xi,X2 =Q,Xklm)\xiX-i-i2X222\ + (1^2, Tj ^ -rj) 



2ny/l - 2773 



exp 



2(1 - 272/3) 



72 + -(2 - 37 )atan 



2-37-2 



4 Unsaddle 



dr] 

The second situation (isotropic Hessian) corresponds to 



(90) 



det V(fe4'" Ud(si)5d(s2) = 



1 xijul - 16{wl+wi)) 
4 (Ai - A2) 



(5d(Ai - A2)<5D(a;2) , 



therefore 
drj 



J dXiPo{'n,\i,X2 = \i) X J dxid'^XkimPi{xi,X2 = 0,Xkim)\xi{ul/4- 4(wl + W2))\ + (1^2,77^-77), 

(91) 



i?2 
1 



TVR 



exp 



n/2^ 



(i + r 



^ Gi{l)Pirl)^ 



7271 i?2 



Both QbIj) and Gsil) a-i'e weak functions of 7 of order unity. The main 7 dependence oc 7"^ reflects R as the fundamental 
scale for the singular points. 

We note that the number density of the "sloping plateaux" is proportional to the density of the saddle points, hence 
this type of singular points is predominantly concentrated near mean fleld values (small 77). In contrast, the number density 
of "bifurcation" points is proportional just to the PDF of the field and, hence, the bifurcation points are as frequent in the 
regions of high field values as in the low ones. This may provide explanation for the observed insensitivity of the curvature of 
the skeleton to the threshold, if we conjecture that most of the curvature accumulates near the "bifurcation" points. 
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2D 3D 



Skeleton 

Anti-skclcton 


4.21R, 
4.21K, 


Skeleton: 

Anti-skeleton: 

Inter-skeleton: 


{4.65R.)2 

{4.65i?*)^ 
(3.38R.)2 


All primary 


2.11R, 


All primary 


(2.36R«)^ 


All secondary 


5.54R, 


All secondary 


(3.02i?.)2 


Total 


1.55iJ* 


Total 


(1.86/?,)^ 



Table 3. Inverse average integrated flux (the characteristic area (3D) or length (2D) per critical line) of the critical lines of different 
types. 



6 CONCLUSION & PERSPECTIVES 



The filamentary structure is a dramatic feature of the observed or simulated Cosmic Web. This paper investigated how 
the set of critical lines of a given field corresponds to an intermediate representation of the field, which is more extended 
than the knowledge of the critical points, but nevertheless much more compact than the field itself. It introduced the stiff 
approximation, which states that the tangent vector to the critical lines only involves up to the the second derivative of the 
fields. Within its framework it has been demonstrated that, for stationary Gaussian random fields, ergodicity allows one to 
recast the description of the ND critical fines into a point process, which only involve the first spectral parameter, 7, when 
considering the differential length as a function of the contrast, and the second spectral parameter, 7, when considering it as 
a function of the modulus of the gradient. The former probability distribution was shown to involve the average flux of the 
Gaussian curvature of the ID sections. In turn, these averages can be carried out analytically almost to the last integral in 
2D and 3D, and provide simple asymptotics at large and small contrast. The detailed contribution of all types of critical lines 



as a function of thresholding was described. The main results of this investigation corresponds to equations (27 1 and (28 1 



for the differential length of the skeleton and the total set of critical lines in 2D and equations (48 1 and (51 1 in 3D. Their 
generalization to N dimensions is given by equation (AI81 in Appendix |A] Table |3] summarizes the average integrated fluxes 
(i.e length per unit volume) of the critical lines. For instance in 3D one expect on average one skeleton line crossing a random 
~ (5-R,)^ surface element. 

These findings were illustrated on scale free power spectra with spectral parameters which are relevant to cosmologjrj 



The prediction of the stiff approximation was checked against measurements for global skeletons ( Sousbie et al. 2008 1 on 
realizations of these fields in two and three dimensions and was found to be in good qualitative agreement. The differential 
curvature of the corresponding lines was also measured (section 5.21 and the corresponding radii were found to be ~ 8-R, and 



~ 2.5-R* near = in two and three dimensions respectively. Hence an access to both the curvature and the length of the 
skeleton provides the means of constraining two shape parameters, 7 and 7. The stiff approximation is also implemented to 



compute the differential curvature in 2D. Finally (section 5.3 1, the stiff theory of the singular points of the critical lines was 



laid out in general, identifying genetically three types of points: critical points of the underlying field, bifurcation points and 
slopping plateaux. Again, the stiff approximation provide means of computing the number density of these points. Appendix [P] 
derived the general joint probability of the field and its successive derivative in arbitrary dimensions, which come into play 
when computing these higher order statistics. 

Clearly the formalism developed in this paper will be useful in the context of the upcoming surveys such as the LSST, or 
the SDSS-3 BAO surveys since it yields access to the shape of the power-spectrum without artifacts related to varying light 



to mass ratio. For instance, Sousbie et al. (2008) first applied the corresponding theory to the SDSS-DR4 catalogue in order 
to constraint the global dark matter content of the universe, since the cosmological parameters are directly a function of the 
spectral parameter, 7. Its application to CMB related full sky data, such as WMAP or Planck should provide insight into, 
e.g. the level of non-Gaussianity in these maps (see SPCNP for a discussion). Similarly, upcoming large scale weak lensing 
surveys could be analyzed in terms of these tools ( Pichon et al.|2009 1 . 

A natural extension of the theoretical component of this work would be to investigate the properties of the bifurcation 
points in anisotropic settings and extend beyond the stiff approximation the preliminary results of Section [5. 3| This will be 
the topic of a forthcoming paper. Another natural venue would be to also investigate the statistical properties of, e.g. the 
peak patch walls (surface, curvature) defined as 3:3 = in the eigenframe of the Hessian. Eventually, a global theory of the 
critical manifolds beyond the local approximation should also be developed to provide a framework to study the connectivity 
of the critical lines. 
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APPENDIX A: THE STIFF SKELETONS OF ND FIELDS 

The emphasis in this paper is on developing the analytical theory of the critical lines of a given GRF in two and three 
dimensions. Yet the critical lines in higher dimensions are of interest in more abstract spaces such as space-time or space- 
smoothing etc. . In 3-1-1 Dimensions, corresponding to 3D space-l-time, the 4D critical lines are the dynamical tracks of 
critical points in 3D. An alternative view is to think of the 4D skeleton as event lines of over densities, while the critical 



points correspond to the position and time of merging events. In fact Hanami (2001) explored sloping saddles (i.e. points in 
position-smoothing space corresponding degenerate saddle points) as a mean of identifying merging events, and argued that 
the ridges (the path of the maxima in position-smoothing space as a function of smoothing) form a 4D skeleton. Clearly these 
higher dimensional spaces would typically not be strictly isotropic, stationary nor Gaussian. As a first step, let us nonetheless 
investigate these N dimensional lines. 

Al Critical Lines in ND 

The local critical lines in A'' dimensions are defined as the points where the condition 

H-Vp = A,Vp, (Al) 
is satisfied. This can be expressed as the condition of the vanishing of the N — 2 antisymmetric tensor 

S = S'i''2>-.'iv-2 = J2 VkpH^e'^ '"-^•''"'VmP = (A2) 

klm 

as defined in equation ([T]). In spaces of dimension > 4 it is more compact to consider the Hodge-dual rank 2 tensor 

(*S) = i*Sf = --^ J2 5'"'"-''"-^^n,...,»«-2.«,. =0 . (A3) 



(N~2y. 
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Local direction of the filament corresponds to the right null-vector Sr*" of the A'^ — 2 + 1-rank tensor of the derivatives of S 



5r = 



(A4) 



A non-trivial solution of this set of C% homogeneous equations generally exists, since the existence of the left null-vector 
X^ii ^iiP (Vfe5*^''^' ' ''""^) = imposes C%_i linear relations leaving exactly C% — C%_i = N — 1 independent equations 
to define a line. 

The notion of primary skeleton lines is automatically generalized for A'^ D as the subset of critical lines obeying 



H-Vp = \i\7p, and Ai + A2 ^ , 



(A5) 



where Ai is the largest and A2 is the second largest of the sorted eigenvalues. 

Let us derive the general expression for the statistical average of the flux of the lines arbitrarily defined over the properties 
of the ND random field by A' — 1 equations — 0, i — 1 ■ ■ ■ N ~ 1, where S^{x, Xk,Xki, ■ ■ ■) are functions of the field and it's 
derivatives. We can shortcut the procedure of flux evaluation by marking each line with one intersection point with a fiducial 
surface E, orthogonal to it, and flnding the N-1 number density of the intersection points on the surface E = 0. The average 
number density of the points defined as the intersection of n non degenerate hypersurfaces a' , . . . , is given by 



J dxd^Xk ■ ■ ■ P{x, Xk,--- )So{(r^) ■ ■ ■ (5D(cr")ldet(Va^, ■ ■ ■ Vct")] 



(A6) 



Let us choose ■ ■ ■ 5*^ ^ as ct^ • • • ct" and E to be . Expanding the determinant |det(VS'^, • • • , \7S^ ^ , VE)j along its last 
row we obtain 



where 



M" = (-1)'=+Met ViS'' 



J dxdJ^Xk ■ ■ ■ P{x,Xk, ■ ■ ■ )<5d(S') • • • 5d(5''^"')<5d(E) 



i = l,...,JV-l 



i = l,...,fc-l,fc + l,...,JV 



(A7) 



(A8) 



are the corresponding minors. By design the E = surface is to be orthogonal to the line and therefore its normal VE and 
the "vector" M = (M*) are parallel. 



= IMIIVEI 



Without loss of generality, we can consider the intersection point to be at r = and take E = ge ■ r where ge is the unit 
vector in the local direction of the filament, hence |VE| = 1. The average (N-l)D number density of intersection points on E 
surface that gives us the average flux C is obtained by integrating the volume number density over the coordinate along ge, 
2 = es • r, with 5d{^) in equation (A7l properly counting exactly one intersection per line 



C = jnd(fi^-r) = j dxd^Xk---P{x,Xk,---)5u{S'-)---Sj,{S^-^)\M\ 



(A9) 



To apply this general formula to the critical lines one must choose an arbitrary subset of TV — 1 linearly independent 



Note that one can also think of C as the average length of lines per unit volume, which is the interpretation we focus on 
in the main text. 



A2 Stiff critical lines in ND 



In the theory of ND critical lines, the N-1 independent functions 5" that define the critical condition ( A2 1 acquire the following 
simple form in the eigenframe of the Hessian of the field 



S — XaX' 



(A, - AO = 0, i / a . (AlO) 



Here a is the index of the Hessian eigenvector that the gradient is aligned with, as is obvious from the solution a;^ = 0, i ^ a. 

In the stiff approximation, the gradients s'^ = Vs' have just two non-zero components, s\ = XiXa{Xa — Ai) (which 
vanishes on the critical line) and s% — XaKi^a — Ai). The vector that determines the direction of the critical line becomes 

^X^-^Xkl[X^l[{\a-X^) ■ (All) 

On the critical line, it has just one non- vanishing component 

M'' =x^-^Y[K{Xa~\i)^\M\ , (A12) 
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dL/d;//P(77)/L 




Figure Al. dC/dr]/ P{rj)/ L in 2D, 3D and 4D as labeled for the spectral parameter values 7 = 0.57,0.65,0.70 {from bottom to top). 
These quantities are derived here by direct numerical integration of equation ( |A18| . The different bundles corresponding to different 
dimensions have been shifted down (by 0.3 for 2D and 0.15 for 3D) for clarity. Note the change in the power of the asymptotic curves. 



which shows that in the s tiff a pproximation we equate the direction of the line with the gradient of the field. Substituting this 
expression into equation (A9l and integrating over = Su{xi) / {xa{Xa — Ai)) we obtain a simple expression for the flux 

of the critical lines in the stifle approximation 



dxdxkiP{x, 0, xu) 



(A13) 



i.e the flux of critical lines (or the length per unit volume) is given by the average absolute value of the Gaussian curvature 
of the field in the space orthogonal to the skeleton. 

Let us write the probability of measuring the set {A^} as 



n d\ Y[{\ - A,) exp (-\q-,{ii, {A.}) 



where Q-y is a quadratic form in Ai and 77 which functional form is 



Q-,{n, {K}) = rf+ ^^l^'^^l^^ + QN{{\i}) .. 



(A14) 



(A15) 



and ni<j(Aj — Aj) is the Jacobian of the transformation to the Hessian eigenframe. Here Qm involves polynomial combinations 
of the eigenvalues of the traceless part of the Hessian (see Appendix [D|: 



QN{{Xij}) = 



N{N + 2) 



'ijXij , With ^ij — ^ij 



(AI6) 



which can be rearranged explicitly in terms of As as: 



QN{{X^}) = {N + 2) 



^(iV-l)^A?-^A.A, 



It now follows that the differential length of the ND-critical lines is for the stiff approximation: 

AT-l 



1 



V ^ I •' •' i^n i<j i>l 



cxp 



(A17) 



(A18) 



Equation (AI81 is the formal generalization of equations (30 1 and (491. For the ND-skeleton, equation (AI8I also holds but 
the integration region should be restricted to the corresponding condition on the sign of the eigenvalues. Since the argument 
of Q-y is extremal as a function of when -yri ~ Ai, the largest contribution at large 777 in the integral should arise when 
Ai oc 777 since near the maximum at high contrast all eigen values are equal ( Pichon fc Bernardeau||l999 1. Hence given that 



ni<j('^« — Aj) is the measure, the only remaining contribution in the integrand comes from |ni>iAi| oc (Ary) , and the 
dominant term at large is given by 



cxp 



1 2 



where Ro — Rt/j is defined in equation ([3|. 
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APPENDIX B: SECONDARY CRITICAL LINES IN 2D 



In this Appendix we present a study of asymptotic behaviour of the lengths statistics of secondary critical lines for 2D 
Gaussian field. Secondary critical lines are the ones that have a gradient of the field aligned with the Hessian eigenvector that 
corresponds to the largest by magnitude eigenvalue, i.e with the direction of maximum curvature of the field. In 2D, this is 
the direction of A2 in the skeleton region, |Ai| < jA2|, and is the direction of Ai in the anti-skeleton region. We shall explicitly 
consider the first type, realizing that the second type is a mirror case with rj —> —r]. 

Our starting point is the part of expression ( |27[ ) that corresponds to the lines where the gradient is aligned with the 
second eigen-direction, in the region when they are secondary, it > 



4^2 



(27r)3/Vl-7^ 



exp \—rf/2\ I dww / du \ 2w — u\ exp 



2(1-72) 



Aw' 



(Bl) 



The absolute value of the transverse to the gradient curvature 2Ai = 2vj — u is evaluated differently for u ^ 2w and u > 2w. 
It is convenient to make the inner integration to be over w, since it can be carried out analytically. The integral splits into 
two terms 



AV2 



drj 



where 
h = 



I2 = 
so that finally 



/ du exp 
Jo 

I du exp 
Jo 



drj 









[ 2(1 




I 

J u 


[ {u- 




/ 


[ 2(1 




Jo 




1 





w{2'w — u)dwe 

S/2 
r-ij/2 

■w{u — 2w)dwe~ 



exp [-V^/2] 



(27r)»/Vl-7 



exp[-r,V2] (/1 + /2) 



16 
16 



/ du exp 
Jo 

oc 

du exp 



(m - 7??)^ 



2(1-7^). 
{u - 77?)^ 



2(1-7^) 



Erfc({i) , 

2 , 

— u- 



Erf (w) 



du exp 



(u - yqf 
2(1-7^) 



|Erf(xi) + J|Erfc({i) 



(B2) 

(B3) 
(B4) 

(B5) 



The integrated length of critical lines U^'^'^ is obtained by marginalization over all threshold values t). Performing this 
integration first 



2^271 Jo 



dfiexp 



u 



^Erf({i) + JjErfc(fi) 



y/2 - ArcCot(2^) 



47r 



= 0.08550 



(B6) 



Thus secondary critical lines are on average almost three times rarer that the primary ones. 



Bl Special cases: rj ^ 00 



At high density thresholds the leading asymptotic behaviour for 70 ^ 1 is obtained by using , ^ = exp — fer^ 

\/2tv{1 — t^) L 



Sd{u — 7?)). Therefore 



dC-^^^o. 1 ^^^pT 2/2. 1 / 2 _ ^ ^g^^ 



B2 Special cases: rj 

At small threshold rj series representation 

(u - jtjf 



exp 



= exp 



2(1-72) 



n!(l 



^ n!(l-72) 



2(1-72) 

where Hermite polynomials H2n are taken in probabilistic notation, gives 
dC 

— 77^ exp [-r;72J 

^ n=0 
1 



2\n/2 



Hn 



(B8) 



977 



^ 00 
(r? 0) = --j= exp [—if /2] ^ An{')ri)" where 



2v/27r7i!(l - 72)"/2 



^ duexp [-u/2] H„(u) ^^/l^u - ^Erf [/T^^u] + ^Erfc [/T^^u]^ 



B9) 
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dL/d7//P(;7) 




-4-2 2 4 



Figure Bl. Differential length of the secondary critical lines dC/dri/PDF in 2D for the complete set (solid) and the ones with Vp 
aligned with A2 direction in the skeleton |Ai| ^ IA2I region (dashed). Different curves from purple to green correspond to the spectral 
parameter values 7 = 0.3, 0.6, 0.95. 



The first three coefficients are 

x/2(l-7') + acot[V2(l-72)] - atan[V2(l - 72)] 



Ao = 



Ai = -r^ 1 + 



V \/2(l-7^) \/(3-272) 

If we add all secondary critical lines, the odd power terms of the expansion ( |B8[ ) cancel, while the even double recovering 
symmetrical behaviour of the differential length with the threshold. In Figure [Bl| this behaviour is illustrated. 

Under our definition of the secondary critical lines, for 7 > l/\/2 there is an excess of critical lines near zero threshold. 
The curvature at 77 = is positive and diverges in the limit 7^1 when our series expansion formally fails. This divergence 
in the second derivative of the differential length is exactly opposite the one the primary lines demonstrate in this limit. We 
should emphasize, that near r; = the behaviour of critical lines of individual type depend significantly on how exactly they 
are defined. 



APPENDIX C: ASYMPTOTIC BEHAVIOUR OF CRITICAL LINES IN 3D 

There are four regions with the different signs of sorted eigenvalues in 3D: I — (0 > Ai ^ A2 ^ A3), II — (Ai ^ 0, > A2 A3), 
III — (Ai ^ A2 ^ 0, > A3) and IV — (Ai ^ A2 A3 ^ 0). Since w is non-negative, the correspondent zones of integration for 
equations (481 and (49l are easy to visualize in (v^u) plane (see Figure Cll. The integration limits and the integrand acquire 
the following form 



/ 
II 
III 
IV 



dw 



dw 



dw 



dv 
dv 
dv 



v-\-Zw 
— v-\-Zw 



2v 
2v 



du 
du 
du 



v — Zw 
v — Zw 



roo rw p — v — Zw 

I dw I dv I du 

•JO -J —w ^ — 00 

Changing variables u - 

dw I dv 



1-2 1-2 

-It V — u 

3 3 

up — -{u + v)'^ + -w{u — 2v) 

w^ — ^({t + v)'^ — ^to(m — 2v) 

1-2 1-2 -2 

— U V — W 

3 3 



- / - 2 -2 \ 

} X w (w ^ V ) exp 



{it-'yrjf 

L '2(1-72) 



V in III and IV one can combine the last two cases with the first two 



I + IV 
II + III 



dw 



dv 



roo 

du 


"1-2 


Sw — v 




pZw — v 




du 


-2 

W — 


2v 





1 -2 -2 



^(u + v)'^ + '^w{u — 2v) 



} X w [w 



-2\ 

v ) exp 



15 _2 

--W 



15 .2 



5-2 



5-2 



(Cl) 

^jiu,v) ■ 
(C2) 
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Figure CI. Left: Integration zones in v, u plane based on the signs of the eigenvalues. Variables are given in units of w. Here v varies from 
—w to +'w, while ii is unrestricted. Three inclined lines are (from top to bottom) a) Ai = => u = —v + SsD, b) A2 = => u = 25 and c) 
A3 = M = — C — 3to. In the upper sector I (that stretches to infinity in ii) (0 > Ai A2 ^ A3), next is the zone II (Ai 0, > A2 A3), 
then III (Ai A2 0, > A3), and, finally extending to minus infinity in u is the sector IV where (Ai ^ A2 5^ A3 0). Centre: two zones 
of integration after variable change leading to equation \C2\ . Horizontal dashed lines mark the further subdivision of the integration 
space if the order of integration is changed according to equation ( |C5| |. Right: Integration zones in the (v,w) plane after u has been 
mapped to the [0 — 00] interval. Variables are given in units of it. The lower triangular zone corresponds to semi-open upper band in 5 — it 
of the centre panel. In this region, the integrand is given by the terms I+III of equation \C2\ . Vice- versa, the open upper band in the 
{v, w) plane corresponds to the II-I-IV integration over the lower triangular zone in the v — u space. The right-most sector of this region, 
however, corresponds to negative ii, so the integrand in this sector has coordinate change u — > — ii,ii — ► —v The dashed lines show the 
subdivided integrals given in equation JCSl, which corresponds to subdivisions in the centre panel. 



where 



$T.(u, rf) — exp 

Direct evaluation of the integrated length gives 



2(1-72) 



+ exp 



(a + 777)2 

2(1-7^) 



L = 0.289627 {xR~' 



(C3) 



To study high threshold regime it is advantageous to make the u integration the outmost one, since it depends on the 
variable threshold 



I + IV : 



II + III 



dw 



dw 



dv du - 
dv j du - 



du I dw I dv + 
Jo J -w 

) fu/2 pSw — u 

du / dw / dv + 

Ju/4 J -w 



i/2 



'L- - 



du / dw 

Ju/4 

ru/2 

dw / dv + 

Ju/2 J -w 



dv 



(C4) 



du 



du 



dw j dv(u,v ^ —u, —v) 

fl/2 Ju/2 



The parenthesis in the last term indicate the substitution that must be performed in the integrand. Right panel in Figure [CT] 
illustrates the integration zones now in the (v, w) plane. 

Although one can perform the v integral analytically and reduce the problem to two-dimensional integration, the resulting 
expression is too cumbersome. We can obtain useful limits already from unreduced formulae. In particular, at high density 
threshold, 77; 00, only the first integral in the term ( |C4[ |, which contains w,v r-^ Q neighbourhood, is not exponentially 
small. Moreover, in the leading order the upper limit of the integral over w can be set to infinity. 



dC 



2r5/2 



3^5 



47rV27r(l - 72) 
1 



exp 



1 2 



oc poo 

du 



POO r 

I dw J 



dvw[w^ — v"^) fit^ — v"^ ~ 3to^1 exp 



{u - jrif 15 .2 5 .2 



27r 



exp 



1 2 



27r 



(C5) 



APPENDIX D: JOINT DISTRIBUTION OF THE FIELD AND ITS DERIVATIVES FOR A GRF 



The joint point distribution functions that are needed for the study of the critical lines in this paper are Po{x,Xki) and 
Pi{xi,Xijk), taking into account that for Gaussian random field there is no cross-correlation between odd order derivatives 
and the field itself or even order derivatives. When considering the curvature of the critical lines, fourth order derivatives, 
and, thus, more general Po{x, Xki, Xkimn) have to be considered. Some well known results in 2D and 3D are first summarized 
in section [Dlj More general results can be obtained by resorting to a general framework which is sketched in section [02] and 
applied in section [D3] for the various cases of interest. 
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Dl Lower order joint distributions 

Distribution of the Gaussian field and its second derivative in 3D. The full expression for Po{x,Xki) for the Gaussian 



field is given in Bardeen et al. ( 19861. Introducing the variables 



u = -Ax = -(xii + X22 + xaa) , w = ^{xii - X33) , u = ^(2x22 - a;ii - 2:33) , (Dl) 

in place of diagonal elements of the Hessian (a;ii, a::22, 2:33) one finds that u,v,w,xi2,xi3,X23 axe uncorrelated. Importantly, 
the field, x is only correlated with it — Ax and 

{xu)='Y, {xv)=0, {xw)=0, {xxki) ^ 0, k ^ I, (D2) 

where 7 is the same quantity as in equation ([5|. The full expression of Po{x, x^i) is then 

Po{x,Xki)dx(fxki = (^2ny/^(l — ^2)1/2 ™P ( ^2 Q2{v,w, a;i2, 2:13, 2:23)] J dx du dv dm dxi2 dxi3 dx23 , 

with the quadratic forms Qo and Q2 given by 

Qo^x^ + ^jT—^ Q2 = 5v^ + 15{w^ + XI2+XI3 + XI3). (D3) 
(1 -72) 

It depends only one a single correlation parameter: 7. 

First and third derivatives of the Gaussian field in 3D. A similar procedure can be performed for the joint probability 
of the first and third derivatives of the fields, Pi(xi, Xijk) by defining the following nine parameters (see also ( Hanami|2001 1): 



Ui = ViU, vi = ~VkVk)x, with j < k, and w^ = ^ l^.Vt ~ '^Aj x , (D4) 

and replacing the variables (aJiii, 0:422, a;i33) with {ui,Vi,Wj). In that case, the only cross-correlations in the vector 
{xi,X2,X3,u-i,vi,wi,U2,V2,W2,U3,V3,W3, X123) which do uot vanish are between the same components of the gradient and the 
gradient of the Laplacian of the field: 

{x,u,) = 7/3, j = l,2,3, (D5) 
where 7 is the same quantity as in equation ([5|. This allows us to write: 

105^'''^3'^ / 1 \ ' 

Pi{xt, Xijk)d'^Xi d^°Xijk = , ^T3/2M -2X3/2 ™P ^ o + Qs) d^Xi d^Ui d^Wi d'^Vi dxi23- (D6) 



(27r)13/2(l_ ^2)3/2 2 

with the quadratic forms: 

\2 



-^J2[ %~J^2) +^?J > Q3 = 105(2.?23 + E("'+»')) • (D7) 



The Gaussian field and its second derivative in 2D. Introducing the variables 

u = -Ax ^ -{xii + X22) , w = ^{xii - X22) , (D8) 
one finds again that u,w, X12 are uncorrelated. The expression for Po(x,Xki) is then 

Po{x,Xki)dxd^Xki = (27r)2(i^ ^2y/2 (^^\ [Qoix,u) + Q2(w,a;i2)]^ dx du dw dxi2 , 
where the quadratic forms Qo and Q2 are 

Qo^x^+ ^"r^f." , Q2 = 8iw^ + xl2). (D9) 
(1-7^) 

First and third derivatives of the Gaussian field in 2D. Defining the following 4 uncorrelated parameters: 

Ui=\/iU, Wi = \/i(\/i\/i-^A) X, (DIO) 



yields 



4 



Pi{xi,Xijk)d^Xid*Xijk = (27r)3^(l^- 7^) '^^^ \ \ + Qs) ) d^Xi d^md^Wi . (Dll) 
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with the quadratic forms: 



Q^-^T.[ %'J^2) ■ Qs = 32Y^wl (D12) 

It is the purpose of the next section to elucidate the nature of these quadratic forms and to show how similar expressions can 
be obtained for any combination of derivatives in a space of any dimension. 



\2 



D2 Theory 

To proceed further, a more systematic way of computing the correlations between the field derivatives is needed. This can 
be provided by the harmonic decomposition of symmetric tensors (such as derivative tensors). The main results are outlined 



hereafter, the reader being referred to Cardoso ( 2009 1 for a detailed exposition 



Harmonic decomposition of symmetric tensors. The harmonic decomposition of symmetric tensors amounts to projection onto 
the irreducible representations of SO(n). It is obtained in close form as follows. A symmetric tensor T of rank n is associated 
with a set {T^^' |0^£^n, n — I even} of "harmonic components" where each T*-^' is a symmetric trace-free tensor of rank 
£. Index £ can be understood as an angular frequency. We refer to it as the "frequency" of the component. The harmonic 
component at frequency ^ = n — 2fc of a rank n tensor is obtained as 

rp(n — 2k) _ rp 

where (tr* • ) means applying k times the trace operator (contraction over any pair of indices) and where T denotes the 
traceless part of tensor T. In indexed notations, the first (ranks 0, . . . , 5) de-traced tensors on are given hy t = t, ti — ti, 

3 6 3 10 5 

with an implicit summation over repeated indices and symmetrization between parenthesized indices (for instance: taa(jSki) = 

[taaj3kl +taak5lj +taalSjk]/S and SO Ou) . 

Invariant statistics. Let T — {To, Ti, . . .} be a set of symmetric tensors which are jointly isotropically distributed. A conse- 
quence of isotropy is frequency decoupling: Ti^' is uncorrelated with T^^ ' if ^ 7^ Further, at any frequency £, the scalar 
product (Ta^^ I Tjj'*') is invariant under rotations. It is convenient to arrange these products at frequency £ into a nii x m£ 
Gram matrix Re where me denotes the number of tensors in T having an harmonic component at frequency £ (this occurs 
whenever rank(r) — ^ is a non-negative even integer): 

where indices a and b run only over the m£ relevant values (the specific ordering does not matter). A further consequence 
of isotropy is that, in the Gaussian case, these matrices form a set of sufficient statistics: the joint distribution of T can be 
expressed as a function of those matrices and nothing else, as seen next. 

Spectral matrices. The 'spectral matrix' Re at frequency £ is defined as the expected value of Re, that is, Re = F,(^Re). For a 
set T of symmetric random tensors with a rotationally invariant joint distribution, one finds 

T'' Cov{Ty^T = J2we tr{ReRi^) , 

e 

where we is a positive scalar, which is equal to 2^ + 1 for tensors in R^. 

Spectral matrices for a GRF. Now, we consider the case when in T — {To, . . . , Tq}, the q-ih tensor Tq is the g-th derivative at 
a given point: ti^^.-.i^ — p/dri-^ ■ ■ ■ dri^ of a stationary random field p with spectrum P(v)- Then T is a set of isotropically 
distributed symmetric tensors and each spectral matrix Re can be expressed as a function of the spectrum. Indeed, if ^ — g 
and t~ q' are non negative even integers, matrix Re has an entry [Re\qq> related to the derivatives of orders q and q' given by 

q — q' 2 
lRe]qq' = (-1) 2 giCTqW > 
2 

with the spectral moments ap defined at eq. Q. The geometric factor ge is the squared ratio ge — (||C^||/||C^||)^ by which the 
norm of the ^-th tensor product of any vector ^ is decreased upon detracing. It is equal to ge = £\/{2£ — 1)!! in dimension 
D = 3. We do not provide explicit expressions for wi and ge in arbitrary dimension since only their ratio we/ ge is needed and 
turns out to have a simpler expression than either we or ge: 

uu ^ {2£ + D-2)U 

ge £\ {D-2)\\ ^'''^> 

Some precomputed values are listed in Table [DT] 
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e = o 


e = 1 


e = 2 


e = 3 


£ = 4 


£ = 5 


D=2 


1 


2 


4 


8 


16 


32 


D=3 


1 


3 


15/2 


35/2 


315/8 


693/8 


D=4 


1 


4 


12 


32 


80 


192 


D=5 


1 


5 


35/2 


105/2 


1155/8 


3003/8 



Table Dl. Values of wt/ge = (2i + D- 2)!!/(£! {D - 2)!!) in dimensions D = 2, 3, 4, 5 for £ 5. 



Summary and rescaled forms. We collect all previous results into a normalized form. Using the normalized spectral shape 

itive it 

X„ = — V"p, 



parameters of def. ([5| and normalized derivative tensors X„ defined as 

1 



dr-i^ ■ ■ ■ dvi. 



one finds that 

X^Cov{X)-'X = Y. ^'^tllo-'m' ^^(^^^^"')' ^'^^ [r,U = <Xf and [r,],, = (-1) ^ 7,,, . (D15) 

Note that the diagonal entries of Ti are always equal to 1. 

Special cases and smaller statistics. Our approach compresses a set of derivative tensors into a set Tt of symmetric matrices 
of size me x me, yielding J]]^ meijni + 1)/2 invariant scalars. There are two special cases where even smaller invariant sufficient 
statistics can be found. 

First, at angular frequency I — the detraced tensors are just scalars so that, for ^ = 0, one has [fo]p, = j Xf'>) = 

Xp^'' x!f\ Therefore To actually is a rank-one matrix: To — vv^ where the entries of vector v are = Xp'^'. Hence, at the 
null frequency, we can further compress the mo(mo + l)/2 statistics (the non-redundant entries of Fq) into mo scalars (the 
entries of v). Of course, the £ = term in the quadratic form also reads: 



tr(roro = v^ro\ 



(DI6) 



Second, there are several cases of interest where me — 2. This happens for instance at ^ = with derivative orders and 2, 
at £ = 1 when considering derivatives of order 1 and 3, at ^ = 2 with derivatives of orders 0,2 and 4, etc. Then, for such an £, 



trff .F7M = tr 



(^a I a) (^a \ 6) 
{b I a) (6 I b) 





■ 1 


-7 


-1 

) 




r'y 


1 





where a and b are rank-^ tensors and 7 is a scalar. Simple algebra yields 

\\b + iaf 



trfF^F,-^ 



\\af + 



1-72 ' 

that is, a form ubiquitous in this paper. However, an equivalent, more regular form is 

1 /„ ,,2 . n,n2^ , 27 



(D17) 



trfFfFrM 



\\af + \\bf) + 



-r 



(a I 6) , 



which has the benefit of stressing that, at such a sufficient statistic is only made of two invariant scalars, namely |ja|j^ -I- ||6||^ 
and (a | 6). In the limit of weak correlation 7 — > 0, one has, of course, tr(f<,F7^) = \\af \\bf . An even more symmetric 
form, which stresses the decorrelation between a -I- 6 and a — 6 is 



tr(ffr,-^) 



||a-h6f ^ ||a- 



■b\\ 



2(1-7) 2(1 + 7)- 



D3 Some applications 

We now work out these expressions in some cases of interest. 



Derivative of orders 0-1-2 in 3D. The case X — {Xo,X2} is the simplest non-trivial case. The theory sketched at sec 
applies straightforwardly. In the not ation s of section D2 , we ar e conc erned with frequencies £ = and 



D2 



1, 102/52 = 15/2 (see table Dll. The quadratic form (D15I then reduces to tr(Foro ^) -I- 
we have here mo 



2 tr(r2r2- 



2 for which, in 3D, 
M. For ^ = 0, 



2 and we can use the specific form (D17l to work out tr(roro ^) with [a, 6] — [Xq"', x'^'^ 



1^, ^aa\ 



that is the ( norm alized) field and the trace of its Hessian. For £ — 2, we have here m2 — 1: we need only scalars. Following 



expressions (D15l again, we have r2 = \\X. 



(2)||2 



= 11X2 



?o + O2 = tr(roro-') + ^ tr(f 2r2-') = 



— XabXab and F2 = (— ) 

{Xaa + ^X) 



(2-2)/2 



2 , 



72,2 = 1. In summary: 
15_ 



(DI8) 
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This is, of course, identical to equation ( D3 1 using the local definitions there. It also shows that the complicated expression 
for Q2 in (D3l is nothing but the the squared Euclidean norm of the detraced Hessian (with a 15/2 prefactor). 

Result for orders 1+3 in 3D. We take X — {XijX^}, that is, the first and third order derivatives of the field. The rescaled 
harmonic components are 



ijk O 



We need frequencies 1=1 and £ = 3 for which, in 3D, w\/g\ — 3, w^/gs = 35/2 (see table Dll. For frequency ^ = 1, we 
have rrii — 2; matrix Fi is 2 x 2 with entries given by equation (D15I, that is, diagonal entries equa l to 1 (as always) and 
off-diagonal entries given by (— l)'^~'^'''^7i,3 = —7. Since Fi is 2 x 2, we can still use equation (D17I and finally obtain In 
summary: 



Wl 



tr(fiFr^) +^tr(f3F3^) = 3 tr 

gi 93 ' 





■ 1 -7- 


-1 - 


{ 


.-7 1 . 





= 3 XiXi + 



XiXiijb XiccXidd 





35 _ _ 


+ 


— XijkXij 


35 




y 


XijkXijk ■ 



(D19) 

-j^ _ ^2 J ' 2 ■"'^''""'J'' ■ (D20) 

This is consistent with equation (D6I and reveals the meaning of 2:123 + +™?) equal to ^XijkXijk i.e. the squared 

norm of the detraced third derivative tensor (with a prefactor 1/6). 

The results for other combinations of derivatives can be derived in the same way. A few results are listed below without 
going into much detail. 



Result for orders 0+2+4 in 3D. We consider X = {Xo,X2,X4}. Hoping to improve clarity, we denote yij — [a^ ^Xj^'jij, 
is, the de-traced contraction of the 4th-order derivative tensor. Explicitly, in 3D: 

1 



that 



yij — Xijaa 



, X aabb^ii 



With this notation and recalling that Xijki denotes the traceless part of Xiji^i (the rescaled 4th-order derivative tensor) 
computed according to the prescription (D13l, the quadratic form is 



X 


t 


" 1 


-7 


7 ' 


-1 


X 


15 


Xaa 




-7 


1 


-7 




Xaa 




Xaabb 






-7 


1 




Xaabb 





{ 


■ 1 


-7 


-1 ■ 




.^7 


1 





Xij yij yij yij 



315 

H — ^ XijkiXijki , (D21) 
o 



where yet another spectral shape parameter has to be defined: 



RR 



(To(T4 RqR* 



77^ 



Needless to say that expression (DI8I obtained for X = {Xo,X2} is recovered by cutting the irrelevant terms from equation 
( D2l| ). 

Result for orders 1+3+5 in 3D. To simplify the notations, we introduce local definitions for the derivative tensors and their 
contractions: 



Va — Xabb , 

and, proceeding as above, we obtain the quadratic form: 

1 -71,3 71,5 "' ^ 
3 tr.; -71,3 1 -73,5 

71,5 -73,5 1 



XaXa Xal/a XaZa 
ya X a ya ya ya Za 
ZaXa Zaya ZaZa 



Xabbc 



35 



tabc Xabcdd j 



1 -73,5 
-73,5 1 



XijkXijk Xijf-iijk 
iij kXijk iij k t-ij k 



693 _ 

H Q XijklmXijklm 
O 

(D22) 



The 2D case. The theory applies to isotropic fields in any dimension. We have already provided expressions for the spectral 
moments Q and the coefficients wg/gi of equation (D14I. It remains to find detracing coefficients. In the 2D case, the first 
(ranks 0, . . . , 5) de-traced tensors on R^ are given hy y = y, y^ — yi, y^j — yij — \yaa5ij, 

3 1 5 5 

Vijk = yi]k-^yaa(iSjk), Vijkl = yt]kl- yaa{ijSkl) + -^ VaabbS {ijSkl) , Vijklm = Vijklm- ^ yaa(ijkSlm) + Y^yaabb(iS]k5lm) ■ (D23) 

For the correlation between the field and its Hessian, we proceed as above in 3D with wo/go — 1 and W2/g2 ~ 4 given in 
table [DT] Therefore the quadratic form is 

Qo + O2 = tr(f oFo"') + 4tr(r2F2"') = x"" + ^"""^^^^^^ + 4 XabXab , 
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in agreement with equation (D9l. For the case of first and third order derivatives, we read wi/ffi — 2 and ws,/gz — 8 from 
table Dl so that, similar to equation (D19l, one finds 



^tr(r.r-) 



H trfrsFg 

93 



2 tr<^ 



= 2 XiXi + 



{ 


■ 1 


-1 ■ 




-7 1 . 





XiXith 



(D24) 
(D25) 



with Xijk ~ X 



■ Xaa(iSjk) SO that SxijkXijk cau be checked to equal Qs in equation (D12l. 



The d-dimensional case. We outhne some results in the d-dimensional case. The de-tracing formulae can be extended to the 
d-dimensional case but, in this paper, we will content ourselves with the correlations between the field and its Hessian: 
X = {Xo,X2}- Therefore, we need only 1 = and £ = 2 so that de-tracing remains trivial: the normalized de-traced Hessian 
given by Xij = xtj — ^ Sij Xaa- Hence for the correlation between the field and its Hessian, we obtain the quadratic form 



X 


t 


■ 1 


~7 


-1 


X 






.~7 


1 




Xaa 



d{d + 2) _ 



(D26) 



which is a straightforward extension of the 3D case of equation ( D18 1 . Just recall that 7 is now defined in terms of the spectral 
moments ([5| and that de-tracing the Hessian requires a factor 1/d instead of 1/3. 
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